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^ ! ABSTRACT 

It is possible that the properties of HII regions during reionization depend sensitively 
on many poorly constrained quantities (the nature of the ionizing sources, the dumpi- 
ness of the gas in the IGM, the degree to which photo-ionizing feedback suppresses 
the abundance of low mass galaxies, etc.), making it extremely difficult to interpret 
\ upcoming observations of this epoch. We demonstrate that the actual situation is 

■ more encouraging, using a suite of radiative transfer simulations, post-processed on 

outputs from a 1024^, 94 Mpc N-body simulation. Analytic prescriptions are used to 

\^ . incorporate small-scale structures that affect reionization, yet remain unresolved in the 

■ N-body simulation. We show that the morphology of the HII regions for reionization 

by POPII-like stars is most dependent on the global ionization fraction Xi. Changing 

O . other parameters by an order of magnitude for fixed Xi often results in similar bubble 

I ■ sizes and shapes. The next most important dependence is on the properties of the 

O \ ionizing sources. The rarer the sources, the larger and more spherical the HII regions 

become. The typical bubble size can vary by as much as a factor of 4 at fixed Xi be- 
^ ■ tween different possible source prescriptions. The final relevant factor is the abundance 

of minihalos or of Lyman-limit systems. These systems suppress the largest bubbles 
from growing, and the magnitude of this suppression depends on the thermal history of 
the gas as well as the rate at which these systems are photo-evaporated. We find that 
neither source suppression owing to photo-heating nor small-scale gas clumping signif- 

■ icantly affect the large-scale structure of the HII regions, with the ionization fraction 

power spectrum at fixed xi differing by less than 20% for /c < 5 Mpc~^ between all the 
source suppression and clumping models we consider. Analytic models of reionization 
are successful at predicting many of the features seen in our simulations. We discuss 
how observations of the 21cm line with MWA and LOFAR can constrain properties of 
reionization, and we study the effect patchy reionization has on the statistics of Lya 
emitting galaxies. 

Key words: cosmology: theory - diffuse radiation - intergalactic medium - large- 
scale structure of universe - galaxies: formation - radio lines: galaxies 
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1 INTRODUCTION 

To interpret existing and future observations of the high 
redshift Universe, we need to understand the morphology 
of the HII regions during reionization. Observations of high 
redshift quasars, gamma ray bursts, and Lya emitters are 
currently probing redshifts at which the Universe may have 



been significantly neutral ([Becker et aP 



I2OO3I : iFan et al.ll2Q06l : iTotani et al.ll2006l : iKashikawal l2006[ 



2OOII: IWhite et all 



ISantos et al. Patchy reionization will leave its 

signatu re in the spect ra of quasars and gamma ray 

bursts (|Haiman Loebl I1999I: iMiralda-Escude et al.1 [2OO0I: 



mmcquinn@cf a. harvard . edu 



iMadau fcR ees 2000l: iFurlanetto et~al]l2004al l and in the cor- 
relation and lumi nosity functions of Lya emitters (|Haimanl 
2002; Santos 20^: IFurlanetto et al.ir2006al l. Starting in 2007, 
the Atacama Cosmology Telescope and the South Pole Tele- 
scope will dissect the high-/ CMB anisotropics. The size dis- 
tribution of HII regions du ring reionization affects the spec- 
trum of these a nisotropics (IMcQ uinn et al ] |2005l : IZahn et al. 
I2OO5I : llliev et al. 2006b ). Finally, 21cm maps of the reion- 
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izing Universe may soon be available. The GMRT, LO- 
FAR, and MWA arrays will begin observing high redshift 
21cm emission within the next few years The 21cm sig- 
nal will be an excellent probe of the structure of reio n- 
ization ([Zaldarriaga et alJ |200 j: iFurlanetto et al.1 l2004bl l3: 



iMellema et al.ll2006l : iFurlanetto et al.l l2006bh 

A proper interpretation of these observations requires 
an understanding of how properties of the ionizing sources, 
how gas clumping, and how source suppression from ther- 
mal feedback impact the size distribution of HII regions. It is 
computationally demanding to simulate reionization in large 
enough volumes to capture the large-scale bubble morphol- 
ogy, and many previous numerical studies simulated only 
a limited number of reionization scenarios, making it diffi- 
cult to isolate the impact of each of the numerous uncertain 
parameters. 

We do not know which objects reionized the Uni- 
verse, although it is most likely that stel lar sources pro- 
duced the bulk of the ionizing photons fe.g. JWvithe Loebl 
(|2003l lV In this case, it is unclear whether the ionizing 
photons were produced by the more numerous galaxies 
with halo masses m ~ 10^ Mq or mainly by rarer, more 
massive galaxies. Locally, the rate at which dwarf galax- 
ies convert gas in to stars scales as g!;ala xy mass to the 
two-thirds power (jKauffmann et all I2OO3I 1. If the same is 
true in the high redshift Universe, then the more massive 
galaxies could dominate the production of ionizing photons. 
However, it might be easier for ionizing photons to escape 
into the inter-g alactic medium (IGM) from smaller galaxies 
(|Wood Loe 3I2OOO). Analytic models predict larger HII 
regions in scenarios in which the most massi ve galaxies pro - 
duce more of the ionizing photons (Furlanetto et al.l 120051 ). 
In spite of our ignorance regarding which sources reionized 
the Universe, numerical studies have yet to examine how 
reionization depends on the properties of ionizing sources. 

Further, we have little observational handle on the 
amount of small-scale structure, or 'gas clumping', in the 
high redshift IGM, and researchers have not reached a con- 
sensus regarding its impact on the morphology of reion- 
ization. Many previous large-scale reionization simulations 
have either entirely ignored structure on scales smaller than 
the simulation grid cell or, despite inadequate resolution, 
have incorporated it via a subgrid clum ping factor calculated 
from their larg!;e v o lume s imulations (Sokasian et al. "2003; 
ICiardi et all |l003|: lliiev'et al. 2006a; Zahn et al. 2006b). 
Recently, there has been some effort to calibrate subgrid 
clumpi ng factors from an ensemble of small-b ox simula- 
tions (iMellema et al.1 12OO6I : iKohler et all l2005l ). However, 
even these efforts are very simplified. No study has tried 
to isolate the effect that gas clumping has on the size dis- 
tribution and morphology of HII regions. If the morphology 
is very sensitive to this clumping, it would be hard to trust 
the results from simulations. 

Another relevant piece of physics is thermal feedback 
from photo- heating the IGM, which can suppress star for- 
mation and potentially alter the morphology of reionization. 
However, the extent to which the structure of reionization 
is affected by such feedback has yet to be adequately ad- 



^ For more information, see |http://www.lofar.org/] and 
|http: / /web, hay stack, mit.edu/ arrays/MWA/[ 



dressed. iKramer et al.l (|2006h , utilizing an analytic model 
for reionization that includes feedback (albeit, on halos that 
cool via molecular line emission), found that it can have a 
dramatic impact on bubble sizes, in some cases creating a 
bimodal bubble size distribution. Similar claims may also 
hold for thermal feedback on galaxies that cool via atomic 
transitions - the more likely culprit to ionize the Universe. 
Ilhev et al ] (2006c) found using radiative transfer simulations 
that thermal feedback plays a key role during reionization, 
marginalizing the contribution from halos with masses below 
10^ Mo. 

In addition, the presence of minihalos and the rate at 
which the gas fron i these halos is p hoto-evaporated may 
shape reionization. Illiev et al. 1 (|2005al ) show that a signifi- 
cant fraction of the ionizing photons will be consumed by 
minihalos and claim that the effect of minihalos on the mor- 
phology of reionization is similar to changing the efficiency 
of the sources. On the other hand, Furlan etto OhI ([20051 ) 
argue analytically that minihalos can create a well defined 
peak in the bubble size distribution that is set by the mean 
free path for an ionizing photon to be absorbed by a mini- 
halo. The effect of minihalos on the characteristics of the 
HII bubbles has not been investigated in simulations. 

In this paper, we present a suite of parameterized mod- 
els, using large volume radiative transfer simulations, to 
understand the impact of each of these uncertain quanti- 
ties on the morphology of reionization. Realistic simulations 
of reionization require extremely large volumes with high 
mass resolution. Previous estimates suggest that, in order 
to capture a representative sample of the Universe during 
reionization, one needs a simulation box with a side length 
of approximately 100 M pc comoving (Barkana Loebll2004l : 
IFurlanetto et al ll2004bl ). To resolve halos at the atomic hy- 
drogen cooling mass (mcooi 10^ Mq at z = 8) in a sim- 
ulation of this volume, one needs about 30 billion particles 
- larger than any N-body simulation to date. In order to 
get around this computational difficulty, we employ a hy- 
brid scheme that combines a 1024*^ particle, 94 Mpc N-body 
simulation with a Press- Schechter merger history tree. The 
merger tree allows us to incorporate halos that are unre- 
solved in our N-body simulation. Additional effects such as 
thermal feedback and minihalo evaporation are incorporated 
in our simulations with analytic prescriptions. 

This paper is organized as follows. In 31] we outline the 
N-body and radiative transfer codes used in this study. The 
radiative transfer code is discussed in more detail in ^A] 
Section ^S) describes our method for including unresolved 
low-mass halos. In 21 we investigate the effect of different 
source prescriptions on reionization, and in JSl we discuss 
the effect of source suppression owing to photo- heating. Sec- 
tion 36] considers the role of quasi-linear gas clumping and 
minihalos in shaping the morphology of reionization. Sec- 
tion 33 discusses the dependence of the morphology on the 
redshift of reionization. The relevance of the previous results 
to observations of Lya emitters and of high redshift 21cm 
emission are discussed in JH] 

Throughout this paper we use a ACDM cosmology with 
Us = 1 erg 0.9, Qrn = Q-3, ^^A = 0.7, Vtb = .04, and 
h — 0.7 (jSpergel et al. | [2003h . All distances in this paper are 
in comoving units. 

More recent measurements suggest that may, in fact, 
be lower than the value assumed in this work ( Spergel et al.l 
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I2OO6I I. The best fit WMAP value is as = 0.74±.05 and when 
combined with other CMB experiments, the 2Df gala xy sur- 
vey a nd the Lya forest becomes ag — 0.78 ± .03 ( Vie l et al.l 
I2OO6I ). A lower as reduces the number of ionizing sources 
during reionization. However, according to analytic models 
for the halo distribution, the sources in a as = 0.8 universe 
are equivalent to those in a as = 0.9 universe at a slightly 
earlier time. Specifically, structure formation in a as = 0.8 
universe at redshift 1 + z should be identical to that in a 
(78 = 0.9 universe at 1 + z' = 9 (1 + z)/S. This occurs be- 
cause halo abundances depend on as through the combi- 
nation D{z) erg, where D{z) ~ 1/(1 + z) is the high redshift 
growth factor. Analytic models for reionization based on the 
excursion set formalism also depend on as only through the 
same combination D[z)aQ. Therefore, if as is lower, this is 
equivalent to a simple re-mapping of redshifts. Furthermore, 
in JTlwe demonstrate that the bubble structure (at fixed ion- 
ized fraction) is relatively independent of redshift and hence 

CTS. 

This paper focuses on predicting the large-scale mor- 
phology of reionization, rather than precisely when reion- 
ization happens. Furthermore, we do not focus on under- 
standing the morphology at times when the global ionized 
fraction is near zero or near unity - in both limits, detailed 
modeling of the complex radiative, thermal and chemical 
feedback processes is essential and challenging. Instead, we 
focus on intermediate ionization fractions. In addition, we 
do not discuss the evolution of the ionizing background or 
the part in lO'^ neutral fraction within the bubbles. We leave 
such discussion to future work. 



2 SIMULATIONS 

We run a 1024*^ N-body simulation in a box of size 
65.6 Mpc with the TreePM code L-Gadget-2 (Springel 
I2OO5I I to model the density field. Outputs are stored on 50 
million year intervals between the redshifts of 6 and 20. 
A Friends- of- Friends algorithm with a linking length of 0.2 
times the mean inter-particle spacing is used to identify viri- 
alized halos. 

The simula t ed h alo mass function matches the 
ISheth TormenI (I2OO2I) mass function for groups with at 
least 64 particles (|Zahn et al ] l2006bl ). However, the mea- 
sured abundance of 32 — 64 particle halos is below the true 
value, but at an acceptable level. Thirty-two particle groups 
correspond to a halo mass of lO^M©. Ideally, we would like 
to resolve halos down to the atomic hydrogen cooling mass, 
^cooi ~ 10^ M0, which corresponds to the minimum mass 
galaxy that can form stars We add unresolved halos into 
the radiative transfer simulation using the prescription de- 
scribed in 33] 

To generate the density grids, we use nearest grid point 
gridding of the N-body particles. Nearest grid point is prob- 
lematic if Poisson fluctuations in the number of particles 
are important at the cell scale. However, a typical cell in 

^ The molecular hydrogen gas cooling channel can lower the min- 
imum galaxy mass. However, Lyman— Werner photons from the 
first stars dissociate the molecular hydrogen, probably eliminat- 
ing this cooling channel prior to the time when the Universe is 
significantly ionized ()Haiman et al.lll997D . 



our fiducial runs has 64 dark matter particles, such that 
Poisson fluctuations are much smaller than the order-unity 
cosmological ones at the cell scale. Nearest grid point affords 
us a higher level of gas clumping (and a more accurate level 
of recombinations) than other gridding procedures, which 
smooth the N-body density field more severely^ 

We use an improved version of the Sokasian et al.l (|200lh 
radiative transfer code, which is discussed in detail in ^Al 
This code is optimized to simulate reionization, making sev- 
eral justified simplifications to drastically speed up the com- 
putation compared to other reionization codes. The code 
inputs the particle locations from the N-body simulation 
as well as a list of the ionizing sources, and it casts rays 
from each source to compute the ionization field. We assume 
that the sources have a soft UV spectrum that scales as zy~^ 
(consistent with a POPII initial mass function (IMF)). The 
parameters we choose for the source luminosities, subgrid 
clumping, and feedback are varied throughout this paper 
and are discussed in subsequent sections. 

The radiative transfer code assumes perfectly sharp HII 
fronts, tracking the front position at subgrid scales This is 
an excellent approximation for sources with a soft spectrum, 
in which the mean free path for ionizing photons is kilopar- 
secs, substantially smaller than the cell size in our radiative 
transfer simulations. 

The radiative transfer simulations in this paper typi- 
cally take two days on a 2.2 GHz AMD Opteron processor 
to reach an ionized fraction of Xi = 0.8. We do not discuss 
ionization fractions larger than Xi = 0.8 in this work because 
our simulation box becomes too small to provide a represen- 
tative picture at larger Xi. In some models for reionization, 
our box is too small even at smaller xi than 0.8 to adequately 
sample the bubble scale and generate clean power spectra. 

We typically choose source parameters so that reion- 
ization ends near z — 7. While overlap - the final stage of 
reionization in which the bubbles merge and fill all space - 
may have occurred at higher redshifts, upcoming observa- 
tions of 21cm emission, QSOs, and Lya emitters are most 
sensitive to low redshifts reionization scenarios. The most 
recent WMAP r = 0.09 =b 0.03 is consistent at the l-a level 
with all the ionization histories in this paper (jSpergel et all 
I2OO6I ). Other papers have attempted to match the source 
propert ies to observations at lower redshifts (e.g., iGnedinI 
(|2000al ll.The escaping UV luminosity of observed galaxies is 
very uncertain, and current observations do not resolve low 
luminosity galaxies at high redshifts. Significant extrapola- 
tion is hence required to connect the properties of observed 
galaxies at lower redshifts to the properties of the galax- 
ies that reionize the Universe. We expect that the source 
prescriptions adopted in this paper are consistent with all 
current observational constraints. 

Table [3] lists the parameters for the reionization simula- 
tions discussed in this paper. A typical luminosity for a halo 
of mass m in the simulations is N(m) = 3x 10^^ m/ (lO^M©) 
ionizing photons . A Salpeter IMF yields approxini ately 
1.5x 10^^ ionizing photons yr (|Hui et al.ll2"ooi ). For 
an escape fraction of /esc = 0.02, for a Salpeter IMF, and for 



^ This is not true for self-shielded regions, which can remain neu- 
tral behind the front (see ^6.2|) . 
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a typical N{m) in our simulations, the star formation rate 
in a halo is S{m) = m/(10^° Mq) Mq yT-\ 



3 UNRESOLVED SOURCES 

Our N-body simulation does not resolve halos with masses 
less than 10^ M©. We use an analytic prescription to include 
smaller mass halos that are sufficiently massive for gas to 
cool by atomic processes and form stars. It is unrealistic to 
ignore the effect of the halos with m < 10^ M©, as many 
previous studies have done, since these halos contain more 
than half of the mass in cooled gas at all relevant redshifts 
(modulo feedback from photo-heating). In addition, halos 
smaller than the cooling mass can still affect the dumpiness 
of the IGM and, thus, are important to incorporate in our 
simulations. 

We outline two methods for adding unresolved halos to 
our simulation in this section and discuss the merits of each 
method. Method 1: We add unresolved halos into each cell 
on the simulation mesh according to the mean abundance 
predicted by Press- Schechter theory. In this text, we use 
this method to include the minihalos. In a cell of mass Mc 
and Unear overdensity today ^o,m, the Press- Schechter mass 
function for halos with mass m < Mc is 



nps(m,5o,M,Mc,2;) 



2 p 



ilogm 



■ So, A 



exp 



2[cr2(m) -cr2(Mc)] 



(1) 



where the function a (M) is the linear-theory variance in 
a region of Lagrangian mass M, p is the mean density of 
the Universe, 6c ^ 1.69/ D(z) , a,nd D(z) is the gijrowth func- 
tion (jPress k Schechter|[T974l : iBond et al.lll99lh . Halos clus- 
ter differently in Eulerian space, and, to account for this, we 
relate the linear overdensity to the Eulerian space over- 
density (5 with the fitting fo rmula calibrated from spherical 
collapse (|Mo White! 1 19961 ): 



So 



1.68647 - 

Sc{z) 
1.68647' 



1.35 



1.12431 



+ 



0.78785 



(1 + ^)2/3 (1 + ^)1/2 (1 + ^) 



(2) 



The radiative transfer code inputs the Eulerian overden- 
sity 6m for all cells from the N-body simulation. To get the 
linear theory overdensity 5o,m we use equation (|2} and Sm- 
In each cell of mass Mc and Unear overdensity ^o,m, we place 
the average number of halos expected from Press-Schechter 
theory using equation ([T]). When including the lower mass 
halos with this method, we need to choose a coarse cell that 
contains more mass than the mass of our largest unresolved 
halo or 10^ M©. We also need a scheme to distribute the 
halos among the cells on the finer grid on which we perform 
the radiative transfer. We discuss this scheme in ^6.21 

The disadvantage of Method 1 is that it involves putting 
the average of the expected number of halos in each coarse 
cell and, hence, ignores Poisson fluctuations in the halo 
abundance. Even the smallest galaxies at these high red- 
shifts are rare and so Poisson fluctuations in their abundance 
can be important. Me thod 2: We account for Poisson fluc- 
tuations by using the ISheth &: LemsonI (| 19991 ) merger tree 




10 100 1000 

M (10 solar masses) 
Figure 1. The halo mass function from our N-body simulation 
and merger history tree algorithm. The merger history tree gen- 
erates all the halos below 10^ Mq, and the more massive halos 
derive from the N-body simulation. We plot curves for the merger 
tree plus resolved N-body halos at z = 6.5 (circles), 8.7 (x^s) and 
11.1 (asterisks). The solid curves are the Press-Schechter, and 
the dashed are the Sheth-Tormen mass functions for these red- 
shifts. The low-mass cutoff for these curves is set by the HI atomic 
cooling mass mcool- The mass function of the merger tree halos 
agree well with the theoretical curves (especially at the lower two 
redshifts, which are the most relevant to this study). 



algorithm to generate the unresolved halos. This algorithm 
partitions a cell with mass Mc into halos and, for a white 
noise power spectrum, produces the correct average abun- 
dance of halos, nps(m,6o,M, Mc, z), as well as the correct 
statistical fluctuations around this mean. The algorithm is 
gua ranteed to work only for a white noise power spectrum, 
but ISheth LemsonI (| 19991 ) show that it works well at re- 
producing nps{m, So,M , Mc, z) and other relevant statistics 
for more general power spectra. This algorithm allows us to 
generate a spatially and temporally consistent merger his- 
tory tree. We find, for the small mass halos of interest, that 
the algorithm generally produces more halos than the Press- 
Schechter prediction. To compensate, we lower as slightly in 
the merger history computation to achieve the best agree- 
ment with the Press-Schechter mass function for our fiducial 
cosmology. 

Figure [1] shows the halo mass function measured from 
our simulations at z = 6.5 (circles), 8.7 (x's) and 11.1 
(asterisks). The merger history tree generates halos below 
10^ Mq, whereas the other, more massive halos are resolved 
in the simulation. The solid curves are Press-Schechter and 
the dashed curves are Sheth-Tormen mass functions for 
these redshifts. The mass function from the merger tree 
agrees best with Press-Schechter and fairly well with Sheth- 
Tormen, particularly at the lower two redshifts - the most 
relevant redshifts for this study. Note that the abundance 
of resolved halos in our simulation is closer to the Sheth- 
Tormen mass function than to the Press-Schechter. 

The merger history tree algorithm generates the halos 
in Lagrangian space, requiring us to then map them to Eu- 
lerian space. The progenitor halos - the halos at the lowest 
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Figure 2. Top Panel: The mass- weighted halo power spectrum 
A^^ for z = 6.6 {thin curves) and 11.1 (thick curves) for several 
source models. [A|^ = (6 ph{k)^) / {27r^), where 5ph{k) is the 
fluctuation in the halo mass density in Fourier space.] The dot- 
dashed curves are for halos with m > 4 x 10-'^'-' M©, the dashed are 
for the resolved halos with m > 2 x 10^ Mq, and the solid are for 
all halos above the cooling mass. We also include the power spec- 
trum of halos above the cooling mass, but weighted as m^/^ rather 
than m {dotted curves) to match a source model discussed in 
All halos with masses between the cooling mass and 10^ Mq are 
incorporated with the merger tree method. Bottom Panel: The 
mass-weighted power spectrum of the source halos (merger tree 
+ resolved halos) at z = 8.7 {solid curve)., and the power spec- 
trum of a smaller box simulation (20 h~^ Mpc, 1024^ particles) 
that resolves halos down to the cooling mass for z = 8.7 {dotted 
curve). The power spectrum from the merger tree method agrees 
with this higher resolution run, boosting our confldence in this 
method. The thin solid curve is k^ h'^pg P^ 5 {k, z) /2tt'^ ^ in which 
is the average Press-Schechter bias and P^s is calculated us- 
ing the Peacock and Dodds fltting formula. 

redshift bin such that they sit on the top of a merger history 
tree - are generated within each coarse cell on a 64"^ grid in 
Lagrangian space, and they are then randomly associated 
with one of the fine cells within its respective coarse cell 
(typically there are 256"^ fine cells). This randomization is 
justified by the fact that Poisson fluctuations dominate over 
cosmological fluctuations at the scale of the coarse cell. To 
map our halos to Eulerian space in a self- consistent man- 
ner, we associate each progenitor halo with a particle whose 
initial (Lagrangian) position is the center of the same fine 
cell as the Lagrangian position of the halo. We then displace 
the particle at each redshift according to second order La- 
grangian perturbation theory (Crocce et al. 2006). At higher 
redshift s, we split the progenitor halo into its daughter halos, 
and all daughter halos are associated with the same particle 
as their parent. This method for adding unresolved halos is 
similar to the PT halo algorithm, an algorithm to quickly 
generate mock galaxy surveys (Scoccimarro & Sheth 2002). 

The bottom panel in Figure [2] plots the mass- weighted 
halo power spectrum A^^ at z = 8.7 from a 1024"^, 20 
Mpc box simulation that resolves halos down to the cooling 
mass {dotted curve). Note that A^ = {5 ph{k)^) / {2ti^) , 



where 6ph{k) is the fluctuation in the halo mass density in 
Fourier space. The A^ of the merger history tree plus re- 
solved halos {solid curve) agrees to better than 40% at all 
scales with the small box A^i^ {dotted curve)^ The level 
of agreement between the dotted and solid curves demon- 
strates that the merger tree method reliably incorporates 
the small mass halos in our simulations. The thin solid curve 
is an analytic prediction for the halo power spectrum given 
by /c"^ P(5(5/(27r^), in which Pss is calculated using the 
Peacock and Dodds fitting formula for the density power 
spectrum, bps — dmmnps{m) bps{m), and bps{m) 

is the Press- Schechter b ias for a halo of mass m [at z = 8.7, 
bps = 3.5] (|Mo &: White! 1 19961 ). This analytic estimate for 
A^^ ignores Poisson fluctuations, and a comparison with the 
other curves indicates that Poisson fluctuations are impor- 
tant on scales oi k > 4h Mpc~^. 

The top panel in Figure [2] shows the mass-weighted 
power spectrum of halos above the cooling threshold from 
the merger history tree method {solid curves) and of the ha- 
los that are well resolved in our box with m > 2 x 10^ Mq 
{dashed curves) at z = 6.6 {thin curves) and z = 11.1 {thick 
curves). The different spectrum of fluctuations between the 
solid and dashed curves suggests that incorporating the un- 
resolved halos may lead to a different HII morphology. As 
the source halos become rarer, their spatial fluctuations in- 
crease and the Poisson component of the fluctuations be- 
comes more important. 



4 SOURCES 

Now that we have a method for incorporating small mass ha- 
los into our simulations, we examine several prescriptions for 
populating the dark matter halos with ionizing sources. We 
consider models where POPII-like sources are responsible 
for the vast majority of the ionizing photons. Even among 
these sources, it is uncertain which galaxies will produce the 
ionizing photons. We consider four models for the source ef- 
ficiencies. In all models, the ionizing luminosity N for a halo 
of mass m is given by the relation N{m) = a{m) m. In sim- 
ulation SI, the factor a is independent of halo mass. Simula- 
tion S2 uses the same source halos as SI except a oc m~^/^ 
(the lowest mass systems are the most efficient at convert- 
ing gas into IGM ionizing photons). In simulation S3, we 
again use the same source halos but set a oc m^^^ (the most 
massive systems are the most efficient). Finally, in simu- 
lation S4, a oc m*^, as in SI, except that only halos with 
m > 4 X 10^° Mq are sources. At 2; = 9, there are 500 
sources in S4 and, at 2; = 7, there are 7000 sources. These 
numbers are in contrast to the other simulations in this sec- 
tion in which there are over 1 million sources at z — 9 and 
over 3 million at 2; = 7. 

Table 1 lists the parameters we use for the runs 
in this section. For simulation SI we set N{m) — 2 x 
10"^^ m/{10^ Mq) photons s"\ To facilitate comparison, we 
normalized the photon production in the S2, S3 and S4 runs 

^ Because the 20 h~^ Mpc box is missing modes larger than 
the size of the box, we expect that it underestimates the true 
spectrum of cosmological fluctuations ( B arkana Sz Loebll2Qo3) . A 
larger box with the same mass resolution would result in better 
agreement between the solid and dotted curves. 
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Table 1. Radiative transfer simulations discussed in this paper. Unless otherwise specified, the subgrid clumping factor Cceii is set to 
unity and the radiative transfer is performed on a 256^ grid. Mq denotes the halo mass in units of 10^ Mq. The functions Cs2, and 

CsA a-re calibrated such that the sources in the respective simulations output the same number of ionizing photons in each time step as 
the sources in simulation SI. 



Simulation Merger Tree Halos* N (photons s ^) Comments 



SI 


yes 


2 X 10^9 Ms 




S2 


yes 


Cs2 M^^^ 




S3 


yes 


CS3 m|/' 




S4 


no 


CS4M8 


includes only m > 4 x 10^° Mq 


Fl 


yes 


2 X 10^9 Ms 


feedback on m < Mj /2; tsf = 100 Myr 


F2 


yes 


2 X 10^9 Ms 


feedback on m < Mj /2; tsf = 20 Myr 


F3 


no 


2 X 10^9 Ms 


includes only halos with m > Mj /2 


CI 


no 


3 X 10^9 Ms 


all cells set to mean density 


C2 


no 


3 X 10^9 Ms 




C3 


no 


3 X 10^9 Ms 


5123 gj.j^ 


C4 


no 


6 X 10'^9 


Cceii given by eqn. |4] 


C5 


no 


6 X 10^9 Ms 


Ccell =4 + 3 Scell 


Ml 


no 


9 X 10^9 




M2 


no 


9 X 10^9 Ms 


Iliev et al. ('2005b) minihalos with mrr^^r^^ > 10^ Mr:^ 

minihalos with mmini > 10^ Mq, fJmh = '^^vir 


M3 


no 


9 X 10'^9 ^g 


Zl 


yes 


1 X 10^0 Ms 




Z3 


yes 







* All radiative transfer simulations are post-processed on a density field that resolves halos down to 10^ Mq. Halo mass resolution is 
extended beyond 109 Mq with a merger tree. Here, 'yes' means the source halo resolution is supplemented with the merger tree down 

to mcool- 



so that the same number of photons are outputted in each 
time step as in SI. In reality, as rarer sources dominate the 
ionizing budget, the rate at which the Universe is ionized 
quickens because the number of high mass halos is grow- 
ing exponentially. Here we are interested in the structure of 
reionization, which is not significantly affected by the dura- 
tion of this epoch. 

The luminosity of our sources only depends on the halo 
mass. This parametrization is most reasonable if, once the 
gas has cooled within a halo, the timescale for its conversion 
into stars is at least comparable to th e duration of reioniza- 
tion ( or a few hundred million years) . ISpringel HernquistI 
(|2003l ) measure a gas-to-star conversion timescale of over a 
gigayear in simulations of high redshift galaxies. However, 
many works in the literature parameterize star formation 
as proport ional to the time deriva tive of the collapse frac- 
tion (e.g., iFurlanetto et al.l (|2004bl lV This parametrization 
assumes that the rate at which a galaxy converts its cold 
gas into stars is much shorter than the duration of reion- 
ization. The effects of alternative parameterizations of star 
formation on reionization are discussed at the end of this 
section. 

The source prescriptions in SI, S2, S3 and S4 are all 
still reasonable in principle. The least massive systems could 
dominate the budget of ionizing photons because it may 
be easier for io nizing photons to es cape from the small- 
est mass halos. IWood Loebl (|2000h find that this is the 
case in static halos owing to the shallower potential well of 
the low mass halos. Internal feedback from galactic winds 
and supernova bubbles may further enhance the escap- 
ing luminosity of smaller halos relative to the more mas- 
sive halos. Internal feedback can also act to shut off star 
formation. ISpringel HernquistI (|2003h find that feedback 



from galactic winds suppresses star formation in the least 
massive systems relative to the more massive. The scal- 
ing a ^ m^/^ taken in model S3 is motivated by the ob- 
ser ved star formation effici ency in low redshift dwarf galax- 
ies feauffman n et al.ll2003l ). 

Because star formation is a complicated process, obser- 
vations rather than theory will likely drive our knowledge of 
the high redshift sources. From present observational con- 
straints, the source prescription used in S4 is closest to being 
ruled out: There is mounting evidence that the highest mass 
halos cannot prod uce enough photons to ionize the Universe 
(|Stark et al.ll2006l ). 

All the simulations in this section were performed on a 
256*^ grid, and the subgrid clumping factor is set to unity 
(i.e., density fluctuations on scales smaller than the cell scale 
are ignored). In subsequent sections, we increase the level of 
clumping and include dense absorbing systems that limit the 
mean free path of photons. Due to the lack of gas clumping 
in the runs in this section, our simulations underestimate 
the number of ionizing photons needed to reionize the IGM. 
However, we find that neither the dense absorbers nor the 
increased clumping have a substantial effect on the topology 
of the bubbles for fixed x^, except in extreme scenarios or 
at higher ionization fractions than we consider. 

Figure [3] compares slices through the ionization field 
from the S2, SI, S3, and S4 simulations (left to right) at 
redshifts 8.7,7.7 and 7.3 (top, middle and bottom panels). 
Panels in a row have the same mass ionized fraction Xi^M- All 
panels have bubbles located around the large-scale overdense 
regions, but the bubbles better trace the overdensities as the 
less massive sources dominate. Reionization in both SI and 
S2 is dominated by the low mass sources and results in a 
nearly identical reionization morphology when comparing at 
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Figure 3. Comparison of four radiative transfer simulations post-processed on the same density field, but using different source pre- 
scriptions parameterized by N(m) = a(m) m. The white regions are ionized and the black are neutral. The left, left-center, right-center 
and right panels are respectively cuts through simulations S2 (a oc m~^/^), SI (o; oc m*-*), S3 (a oc m^/^), and S4 (a oc m*-*, but only 
halos with m > 4 x 10"'^'-' Mq host sources). For the top panels, the volume ionized fraction is Xi^v ~ 0-2 (the mass ionized fraction is 
Xi,M ~ 0-3) and z = 8.7. For the middle panels, Xi^v ~ 0-5 {xi^m ~ 0-6) and z = 7.7, and for the bottom panels, Xi^v ~ 0-7 (^i,M ~ 0-8) 
and z = 7.3. Note that the S4 simulation outputs have the same Xi^M-, but Xi^v that are typically 0.1 smaller than that of other runs. In 
S4 the source fluctuations are nearly Poissonian, resulting in the bubbles being uncorrelated with the density field (xi^y x^^m)- Each 
panel is 94 Mpc wide and would subtend 0.6 degrees on the sky. 



fixed Xi. The HII regions in S3 are larger and more spherical 
than they are in SI and S2. The bubbles are still larger in 
S4. 

The differences between the ionization maps owe to the 
bias differences between the sources. As the sources become 
more biased, they become more clustered around the densest 
regions, resulting in the bubbles becoming larger. In Press- 
Schechter theory, the luminosity-weighted average bias at 
2; = 8 is bps = 2.8 for the S2 source prescription, 3.2 for the 
SI, 5.0 for the S3, and 7.3 for the S4. The S4 sources are 
located in the highest density peaks in the Universe and the 
fluctuations in the density of these sources is the largest. [See 
respectively, the thin solid, thin dotted, and thin dash-dotted 
lines in Fig. [H for a comparison of the luminosity-weighted 
power spectra for the SI, S3, and S4 sources at z = 6.6.] 
The differences between the ionization maps for SI, S3, and 



S4 should allow observations to distinguish between these 
scenarios (as discussed in 31]). 

This trend of bubble size increasing wi th average source 
mass was predicted in the analytic work of iFurlanetto et aP 
(|2005l ). Analytic models typically ignore Poisson fluctua- 
tions in the source abundance, which can dominate over cos- 
mological fluctuations when relatively massive sources dom- 
inate the photon productionlf] For example, the bubble scale 
in S4 is roughly 20 Mpc at 2; = 7.7 and Xi = 0.5 - a scale 
where Poisson fluctuations dominate over the cosmological 
ones. For the SI and S2 source models, cosmological fluc- 
tuations dominate over Poisson fluctuations on the scale of 

^ While it is difficult to incorporate Poisson fluctuations 
in analytic mo dels based on the excursion set formalism, 
IFurlanetto et al . (2005) investigated the effect of Poisson fluctu- 
ations using such models. 
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R (Mpc/h) 

Figure 4. The volume-weighted bubble radius PDF for the SI 

(solid curves), S3 (dot- dashed curves), and S4 (dotted curves) 
simulations. See the text for our definition of the bubble radius 
R. We do not include curves for the S2 simulation because they 
are similar to those for SI. The thin curves are at z = 8.7 and 
Xi^M = 0-3, and the thick curves are z = 7.3 and xi^m = 0-8- 
Simulation S4 has the rarest sources and the largest HII regions 
of the four models. 




0.1 1 10 

k (h Mpc'b 

Figure 5. The ionization fraction power spectrum /S.xx(k)'^ = 
Pxx(k)/27T'^ for the SI (solid curves), S2 (dashed curves), S3 
(dot- dashed curves) and S4 (dotted curves) simulations. For the 
top panels, Xi^v ~ 0-2 (xi^M ~ 0-3), for the middle panels, Xi^y ~ 
0-5 (xi^M ~ 0-6), and for the bottom panels, Xi^y ~ 0-7 (xi^M ~ 
0.8). In all panels, the fluctuations are larger at /c < 1 /iMpc"-*^ in 
S3 and S4 than they are in SI and S2. As the most massive halos 
contribute more of the ionizing photons, the ionization fraction 
fluctuations increase at large scales. 



a typical bubble, but Poisson fluctuations can be important 
in small er bubbles. This de ficiency of analytic models was 
noted in 'Zahn et aD (|2006b"). 

Figure Hlplots the bubble size distribution from 81 (solid 
curves), S3 (dot-dashed curves), and S4 (dotted curves) at 
Xi^M — 0.3 (thin curves) and Xi,M = 0.8 (thick curves). The 
S2 simulation is not included here; it yields bubble sizes that 
are similar to those in SI. How do we define the bubble "ra- 
dius" since the bubbles are far from spherical? For each cell 
in the box, we find the largest sphere centered around this 
cell that is 90% ionized. We say that each cell is in a bubble 
of size equal to the radius of this sphere. We then tabulate 
the radius from all the ionized cells to calculate the volume- 
weighted bubble PDF (zero-radius bubbles are not included 
in the tabulation). This definition of bubble size is chosen 
to facilitate comparison with analytic models of reioniza- 
tion based on the excursion set formalis m in which the bub- 
ble radius is similarly defined (|Furlanetto et al.ll2004bl ) . The 
bubbles are largest in S4 and smallest in SI, and in all runs 
there is a characteristic bubble radius. 

It is useful to compare the measured bubble size distri- 
bution to the size distribution predicted in analytic models. 
The "log-normal" distribution of bubbles found in analytic 
models is present in these simulations. The bubble size dis- 
tribution becomes more sharply peaked in \og(R) with in- 
creasing Xi in our simulatio ns, a trend th at was predicted 
by analytic models (Furlane tto et al.ll2005l V A more detailed 
comparison of the bubble size s between t h ese sim ulations 
and analytic models is given in IZahn et al. 1 (|2006bl ). 

Figure [5] plots the dimensionless ionization fraction 
power spectrum A^^, for the four simulations (SI - solid, 
S2 - dashed, S3 - dot-dashed, S4 - dotted) for the vol- 
ume ionized fractions Xi,v ~ 0.2 (top panel), Xiy ^ 0.5 
(middle panel) and Xiy ^ 0.7 (bottom panel). [The Xi^v 
for S4 is ~ 0.1 smaller than these values.] Note that 
A,,(A:)2 = A:2p..(A;)/(27r2) in which (27r)^ P,,.(A;) fc(k - 
k^) = {xi(\^)xi(\<i)) . For some A^^, the power peaks at the 
box scale (k ~ 0.1 Mpc~^), particularly at larger x^. This 
indicates that there is substantial power in ionization frac- 
tion fluctuations on scales larger than our simulation box in 
some of the considered models. Therefore, the box we use is 
too small to make statistical predictions about reionization 
for some of the models and at some Xi. Lyman- limit systems 
or minihalos may reduce the size of the largest bubbles and 
alleviate this difficulty (see ^6.2|) . 

It is useful to note that an ionization field that is com- 
posed of fully neutral and ionized regions with total ionized 



fraction xiy has variance of Xi^ 
implying that 



on small scales, 



f 

Jo 



dlogk Axx(k)'^ = Xi,v 



-2 
- Xi y. 



(3) 



Because of equation (|3]) and because the snapshot from S4 
has more power on large scales, the snapshots from SI, S2, 
and S3 must have more power than S4 on small scales for 
the same Xi. The distribution of power has important im- 
plications for upcoming observations. Generally speaking, 
the more power on large scales (k < Ih Mpc~^), the more 
observable the signal (see 33) • 

The picture of reionization seen in simulations SI, S2 
and S3 is differen t from that seen in the simulations of 
llUev et al ] (|2006al ). Their simulations resolve halos with 
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m > 2 X 10^ Mq, and reionization ends at z ~ 12 in their 
calculations. Hence, the resolved halos in their simulations 
are very rare and, of the four source models we consider, are 
most similar in abundance to the source halos in S4. Their 
reionization snapshots give the visual impression of many 
overlapping spheres. We do see, particularly in simulation 
S4, that the bubb les become raore sp herical as the sources 
become rarer. See IZahn et al.1 (|2006bl ) for further compari- 
son. 

The prescription we use for the luminosity of the sources 
is simplistic. In all of our source models, the luminosity of 
a halo is monotonic in the halo mass such that the charac- 
teristic source mass is either ttIcooI or m* - the halo mass 
that characterizes the transition to the exponential tail in 
the luminosity function. Star formation is complicated, and 
the characteristic mass of a source could be an intermediate 
mass between mcooi and m*. In this case, the bias of the 
sources will fall between the source bias in S2 and in S3, 
and, therefore, the bubbles sizes will be between the sizes in 
S2 and in S3 if we compare at fixed Xi. 

Surely the luminosity of galaxies depends on additional 
parameters besides the halo mass. Other studies have pa- 
rameterized the luminosity of the sources as proportional 
to the time derivative of the collapse fraction, considering 
the accretion of gas onto sources as a better proxy for the 
star formation rate than the gas mass of the sources. We 
have run simulations with the luminosity proportional to 
the time derivative of the collapse fraction in a cell. We find 
that the morphology of reionization is very similar between 
this parameterization and that of the constant mass-to-light 
model. The reason for this similarity is that the collapse frac- 
tion in a given region is changing nearly exponentially with 
time and so the rate of halo mass growth is proportional to 
the halo mass. Alternatively, star formation or quasar ac- 
tivity may be correlated with major mergers (see Hopkins 
et al. 2006a, b, Li et al. 2006 for discussion). Since major 
mer ger events are more b iased, this results in larger bub- 
bles. Cohn &: Chand (|2006h used an analytic model to derive 



the bubble sizes in merger-driven scenarios. In addition, it 
might have been possible for the gas in smaller mass galax- 
ies (m > 10^ ^o) to cool via H2 transitions. If this is the 
case, stars would form in halos with smaller masses than are 
considered here. These sources would be less biased, and, 
therefore, the HII regions would be smaller and more frag- 
mented Q 



5 SOURCE SUPPRESSION FROM 
PHOTO-HEATING 

The extent to which photo-heating from a passing ioniza- 
tion front affects the ionizing sources and, as a result, the 
reionization process is not well understood. Often, when in- 
cluded in a study, the effect of photo- heating is parameter- 



^ If molecular hydrogen cooling does happen at low redshifts, 
then it may occur in halos with m ~ 10^ ^0- Feedback processes 
may destroy the H2 in smaller halos. However, bps = 2.6 for halos 
with m > 10^ Mq at 2; = 8, as opposed to 6ps = 2.8 in S2, 
such that the bubble sizes will be similar to the sizes in S2. The 
harder spectrum of POPIII stars will make the ionization fronts 
less sharp. 



ized in a simplistic fashion: Star formation is assumed to be 
completely shut off in the low mass sources as soon as an 
ionizing front has passed. However, sources that form prior 
to a front passing will have a cool reservoir of gas with which 
to make stars. Since photo-heating can suppress subsequent 
accretion onto these objects, eventually this reservoir will 
run dry and all the gas will have been converted to stars. 
The timescale over which this reservoir will be depleted is 
uncertain (see discussion in 311) • 

Furthermore, the mass threshold at which sources will 
be suppressed by photo- heating is fairly unconstrained. Of- 
ten the suppression mass scale is taken to be the linear the- 
ory Jeans mass Mj. This choice is, however, problematic. 
The gas will not instantaneously respond to photo-heating 
- there will be some delay, leading to a time dependent sup- 
pression threshold that only asym ptotically approaches the 
Jeans mass for linear fluctuations (|Gnedin HuillToOSh . In 
addition, a spherical perturbation that collapses at z = 8.0 
was at turnaround at z = 13.3. An ionization front passing 
this collapsing mass at, say, 2; = 9, will do little to pre- 
vent the gas from cooling. The collapsing gas is already 
significantly overdense prior to front- crossing, giving it a 
large c ollisional cooling rat e and possibly allowin g it to self 
shield (jPiikstra et al.ll2004l ). lDiikstra et al.1 (|2004h finds in 1- 
D simulations that a substantial fraction of collapsing den- 
sity peaks with mass below the Jeans mass threshold (or 
2.7 X 10^ M© at z = 7 for Tgas = lO'^ K) are still able 
to collapse and f orm g as- ri ch halos in ionized reg ions, and 
iKitavama et al.l (|2000h and Kitavama et al. J200lh find an 
even larger fraction than lPiikstra et al.l (|2004[ ) in 3-D simu- 
lations. 



Illiev et al ] (|2006d ) was the first to investigate with 
large-scale simulations of reionization the effect feedback on 
the sources from photo-heating has on reionization. They 
applied the rather extreme criterion that star formation in 
all halos below 10^ Mq is shut off after 20 million years in 
ionized regions. They concluded from this study that the 
small halos do not play an important r ole in ionizing the 
IGM. Here we expand upon the work of Illiev et al.l (|2006cl ) 
to include more general parameterizations for the feedback 
from photo-heating. 

The parameterizations we adopt for source suppression 
owing to photo-heating are simplistic. However, we show 
that the structure of reionization is largely unaffected by 
feedback even for an aggressive parametrization of suppres- 
sion. If an ionizing front passes a source with luminosity 
Lq at time ti then at time t we set its luminosity to be 
L{t) = Lq exp[— (t — ti)/TsF], where tsf can be thought of 
as the timescale over which the cool gas in the potential well 
of a source is converted into stars. We set tsf — 100, 20 and 
million years in simulations Fl, F2, and F3, respectively. 
We assume that this luminosity suppression affects halos 
with masses below Mj/2, where Mj is calculated for gas at 
T = 10^ K (or 1.4 X 10^ Mo at z = 7). This fixed suppression 
mass misses the time dependent response of the gas to photo- 
heating. The suppression mass Mj/ 2 is approxima t ely th e 
suppression mass found at 2; = 7 in IPiikstra et all (|2004 ) . 
This suppressi on mass is an or d er of magnitude larger than 
that found by IKitavama et al.1 (|2000h . Furthermore, we as- 
sume that halos that form in already ionized regions with 
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Figure 6. The impact of thermal feedback on the ionization his- 
tory. The ionization history Xi^Y(z) for the SI (solid curve), Fl 
(dotted curve), F2 (dashed curve), and F3 (dot- dashed curve) sim- 
ulations. Simulation SI has no feedback, and F3 has maximal 
feedback. Feedback extends the duration of reionization by as 
much as 200 Myr in our simulations. 



SI Fl 




Figure 7. The impact of thermal feedback on the morphology of 
reionization. Slices through snapshots in which Xi^v ~ 0-7. The 
white regions are ionized and the black are neutral. SI includes no 
thermal feedback, Fl has minimal feedback with tsf = 100 Myr, 
F2 has strong feedback with tsf = 20 Myr, and F3 uses only 
halos with m > Mj/2 as the sources (or, equivalently, tsf =0). 
Notice that even the ionization maps from SI and F3 do not 
differ by much, which shows that feedback has a small effect on 
the structure of reionization. 



masses below Mj/2 have zero star formation and do not 
contribute to reionizationQ 

First, in agreement with previous studies, we note that 
thermal feedback can delay and extend the reionization pro- 
cess (Fig. [6|. Simulation SI (solid curve) does not include 
feedback, whereas simulation F3 (dot-dashed) includes max- 
imal feedback (tsf — 0). The duration of reionization is ex- 
tended by about 200 million years in this case. For the other 
feedback scenarios (Fl and F2), reionization is extended by 
a shorter period (100 and 150 million years). 

Figure [71 displays slices through snapshots with = 
0.7 for the simulations SI, Fl, F2, and F3. In SI, halos be- 
low the mcooi always contribute ionizing photons, whereas in 
simulation F3 only halos above Mj/2 contribute. The differ- 
ences between SI and F3 are minor: The small mass sources 
do not change the structure of reionization significantly. SI 
has more small bubbles, and the HIT fronts have more small- 
scale features. Simulation Fl (tsf — 100 Myr) is most sim- 
ilar to SI - long gas-to-star formation timescales essentially 
negate the effect of feedback, and simulation F2 (tsf — 20 
Myr) has less structure in the voids than Fl. In conclusion, 
simulations SI and Fl-3 have a very similar morphology at 
fixed Xi. Feedback does not significantly affect the structure 
of reionization. We find that this conclusion still holds if we 
compare at other Xi^v as well. 

To make the comparison of feedback models more quan- 
titative, we contrast the A^cc at Xi = 0.7 for these four mod- 
els. We find that the A^x of the SI, Fl, and F2 models agree 
to approximately 10% at all scales and that the /\xx of the SI 
and F3 models (no feedback and maximal feedback models) 
differ by at most 20%, with the largest differences being for 
modes near the box scale and for modes with k > 5h Mpc~^ . 

It is simple to understand why thermal feedback has 
little impact on the size distribution and morphology of HII 
regions (provided we compare at fixed Xi). The bubble size 
distribution and morphology are mainly sensitive to the bias 
of the ionizing source host halos and to Poisson fluctuations 
in the halo abundance for sufficiently rare source halos. The 
top panel in Figure [2] compares the luminosity- weighted halo 
power spectrum A^f^ for halos above the cooling mass at 
z = 6.6 (thin, solid curve) compared to A^;^ for halos with 
m > 2 X 10^ Mq (thin, dashed curve). Notice that the differ- 
ence between these curves is less than the difference between, 
for example, these curves and those for the S3 sources (thin, 
dotted curve). In terms of the Press-Schechter bias at ^ = 8, 
^ps = 3.2 for the SI sources whereas bps = 4.3 for halos with 
m > Mj/2 = 1.4 X 10^ Mq. These values should be con- 
trasted with bps = 5.0 for the S3 sources and bps = 7.3 for 
the S4 sources. Therefore, if halos with m < 1.4x10^ Mq are 
evaporated (as in this section), the morphology of reioniza- 
tion is not changed as substantially as the difference between 
the morphology in the SI and in the S3/S4 simulations. In 



^ For simplicity, we take the sources that exist with masses below 
Mj/2 at the instant a region becomes ionized to be the sources 
with m < Mj/2 that contribute photons for all subsequent times. 
In reality, a fraction of these halos that form prior to front crossing 
will become incorporated in more massive halos than Mj/2, halos 
we count as separate sources. Therefore, we double count some of 
the mass in halos and underestimate the effect of feedback when 
TSF > 0- This underestimate does not change our conclusions. 
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fact, Figure [7] shows that the bubbles are largely unchanged 
by feedback. 

All simulations in this section are parameterized such 
that N oc m and such that the suppression scale is Mj/2. 
For lower Xi^v than are shown in Figure [71 the effect of feed- 
back in our simulations is even less significant. If the highest 
mass sources are more efficient at producing ionizing pho- 
tons, reionization will be extended by a smaller amount by 
feedback than we find, whereas if the low mass sources are 
more efficient, feedback will extend reionization by a larger 
amount. The conclusion that the structure of reionization is 
only modestly affected by feedback holds even if the sources 
near mcooi are more efficient at producing ionizing photons 
then we have assumed: We found in 311 that as we made the 
low mass sources more efficient, the properties of the HII re- 
gions are essentially unchanged (compare the panels from SI 
and S2 in Fig. [3]). Lastly, we believe that our choice of Mj/2 
is a fairly extreme suppression mass for low redshift, POPII 
star reionization scenarios owing to effects mentioned at the 
beginning of this section. If the suppression mass is larger 
than Mj/2 or if reionization happens at a higher redshift 
but with the same suppression mass, thermal feedback wih 
be mo re in aportant. However, at z > 10 both lPiikstra et aP 
(|2004l ) and lKitavama et al.1 (|200Q| ) find that the suppression 
mass is much lower than 10^ Mq . 

If molecular hydrogen cooling is able to cool the gas 
in a halo to form a galaxy then most star formation could 
take place in halos with m ^ mcooi- In such a case, thermal 
feedback could play a mo re important r ole in shaping the 
structure of reionization. iKramer et al.l (|2006l l found that 
this scenario could lead to a bimodal bubble size distribu- 
tion. (Note that in the models that we consider in which only 
halos with m > mcooi form stars, feedback does not create 
a bimodal bubble size distribution, and the size distribution 
of the bubbles is largely unchanged by thermal feedback.) 



6 EFFECT OF DENSITY INHOMOGENEITIES 

Density inhomogeneities potentially play an important role 
in shaping the HII regions during reionization. On small 
scales, density inhomogeneities lead t o the outside- in reion- 
ization observed in the simulations of iGnedinI (|2000al V The 
role of these inhomogeneities on the large-scale bubble mor- 
phology has not been investigated in detailed simulations. 
Analytic models make simplistic assumptions to incorporate 
their effects. These models spherically average the density 
fluctuations in a bubble and typically treat a higher level 
of recombinations as equivalent to decreasing the ionizing 
efficiency of the sources. 

Previous large-scale radiative transfer simulations of 
reionization either ignored subgrid density inhomogeneities 
entirely, or they calibrated their subgrid clumping factor 
from smaller simulatio ns (iMellema et al.ll2006l : iKohler et al.l 
I2OO5I ). A simulation of iMellema et al.1 (|2006l l uses a clump- 
ing factor th at is independent of 5 and Xj and neither the 
simulations of 'Melle ma et al.1 (|2006') nor lKohler et al.1 (|2005l ) 
include a dispersion in the clumping for a cell of a given 
overdensity. Both studies of clumping also assume that the 
clumping factor is independent of the local reheating and 
ionization history, which is incorrect in detail. In linear the- 
ory, the smallest gas clump - which is intimately tied to 



the gas clumping fac tor - is given by the filtering mass Mf 
(jOnedin Hui|[l998l V and this mass incorporates the time- 
dependent gas response to heating (see ^C}. The filtering 
mass provides some framework to understand the small-scale 
gas clumping. It is important to understand how sensitive 
the characteristics of reionization are to gas clumping - to 
what extent can gas clumping be ignored or included in only 
a primitive manner? 

Minihalos - virialized objects with Tvir < 10^ K - con- 
tribute to the clumping differently than does the diffuse 
IGM. These virialized objects are unresolved in all current 
large-scale reionization simulations. Minihalos are extremely 
dense and act as opaque absorbers until they are photo- 
evaporated. Since the inner regions of minihalos are self- 
shielded, it is difficult to describe the effect of minihalos 
with a subgrid clumping factor. In addition, most photons 
that pass through a cell should not be affected by a minihalo 
because the mean free path for a ray to intersect a minihalo 
can range between 1 and 100 Mpc. Absorptions by minihalos 
are unimportant when the HII regions are much smaller than 
the mean free path. Once the bubble size becomes compara- 
ble to the mean free path, minihalos may be the dominant 
sinks of ionizing photons within a bubble. IPurlanetto OhI 
predict that minihalos create a sharp large-scale cut- 
off in the size distribution of bubbles, particularly when the 
Universe is largely ionized. If this prediction is true, large 
scale topological features during reionization can be used to 
probe small-scale density fiuctuations. 

We split the discussion in this section into two com- 
ponents: (1) quasi- linear IGM density inhomogeneities, and 
(2) the minihalos. (Our discussion on the effect of minihalos 
also applies to the effect of a more general class of dense 
absorbers, Lyman- limit systems.) The technology needed to 
describe these two forms of gas clumping is quite different. 
In this section, we use only the halos that are well resolved 
in the simulation as our sources (m > 2 X 10^ Mo), and we 
set the luminosity proportional to the halo mass. While this 
source prescription is probably unrealistic, we found in JSl 
that including less massive halos does not change consider- 
ably the structure of reionization. 

6.1 IGM clumping 

We cannot reaUstically calculate the dumpiness of the gas 
from the N-body simulation used in this paper. In order 
to investigate the effect of the clumping, we consider four 
toy models for clumping of the IGM. Simulation CI uses a 
256^ grid, setting the baryonic overdensity to zero and the 
subgrid clumping factor Cceii to unity in every cell. In other 
words, the IGM is completely homogeneous in this model. 
Simulation C2 is a 256^ simulation also with Cceii = 1, but 
it uses the gridded N-body density field. The cell mass in C2 
is 2 X 10^ Mq, approximately the Jeans mass for lO'^ K gas 
at z = 6. Simulation C3 is a 512*^ simulation with Cceii = 1- 
The cell mass in C3 is 3 x 10^ Mq, below the Jeans mass 
at relevant redshifts, but possibly above the filtering mass. 
Table 1 lists the specifications used in the Cl-4 simulations. 

When the Universe becomes reionized, the filtering 
mass Mf can be orders of magnitude smaller than the Jeans 
mass. It takes hundreds of millions of years for the gas to 
respond fully to the photo-heating and clump at the limiting 
scale. Therefore, the 512^ run is closer to reality than the 
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256^ one, but still underestimates the effect of clumping on 
the IGM. To account for this higher degree of clumping, we 
run simulation C4. This is a 256"^ simulation with twice the 
ionizing efficiency of the other runs such that overlap occurs 
at around the same time. In addition, we set the subgrid 
clumping factor in C4 to 
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where kf is the scale that contains the mass Mf (which is 
given by equation IC1|) , Wceii is the cell window function, 
and /cs is the wavevector that corresponds to 10^ M© at the 
mean density - the minimum mass baryonic clump that we 
allow, consistent with a minimal amount of reheating. For 
simplicity, we use a spherical top hat in real space that has 
the same volume as a grid cell for Wceii- We use the Peacock 
and Dodds power spectrum for P55{k,z). The filtering mass 
Mf depends on the redshift at which the cell was ionized. 
Once a region is ionized, this mass increases with time and 
Cceii typically decreases. 

Equation (|4]) would be correct if the window function 
of a cell were instead a top hat in /c-space, if mode coupling 
were absent between modes smaller and larger than the cell 
scale, and if the quantity Mf were appropriate outside of 
linear theory (there is evidence that it is appropriate [^C]). 
Since we are considering non-linear scales, mode coupling is 
important and tends to make the more massive cells have 
higher clumping factors than equation Q predicts. In the 
limit in which most of the density fluctuations are at scales 
smaller than the cell, equation ^ predicts that the number 
of recombinations (oc Cceii Pceii Pceii) is independent of the 
cell's density. This prediction is probably unphysical. 

Note that we assume that the gas clumping in a cell 
is independent of the cell's ionization fraction in all of the 
simulations. This assumption is justified for the gas in the 
diffuse IGM because this low density gas stays almost fully 
ionized when an ionization front passes, provided that there 
is an ionizing background. Virialized objects, such as mini- 
halos, in which the local ionized fraction can be a function 
of density, are included in the computation in ^6.21 

The reionization scenarios in this section reach Xi = 0.5 
near z = 7. The reionization epoch in simulation C4 is 
slightly more extended than the other scenarios owing to an 
enhanced number of recombinations. The volume- averaged 
clumping factor in ionized regions Cv is 30 at ^ = 7 in 
C4, whereas it is 1.6 in C2 and it is 2.7 in C3. The total 
number of IGM photons that escape into the IGM per ion- 
ized baryon is 3 in C4 at the end of reionization, whereas 
it is between 1.2 — 1.3 in C2 and C3. (The recombination 
rate is proportional to the clumping factor.) Note that we 
have removed the particles that are associated with halos 
from the density grid in these simulations since most ab- 
sorptions within these halos are already encapsulated in the 
factor /esc- Other studies left the halos in the density field 
(|Gnedinll2060al : ICiardi et al.ll2003h . yielding a large number 
of recombinations within the source cells (which can be at 
hundreds of times the mean density) and therefore a larger 
photon to ionized baryon ratio. 

Figure [8] depicts a slice through the box at Xi^v ~ 0.5 
for the CI, C2, C3 and C4 simulations. The ionization field 




Figure 8. The impact of gas clumping on the structure of reion- 
ization. A slice through the CI (top left), C2 (top right), C3 (bot- 
tom left) and C4 (bottom right) runs at z ~ 7 and Xi ~ 0.5. All 
the cells in CI are at the the mean density. Simulation C2 is run 
on top of the N-body simulation density field gridded to 256^, 
and C3 is the same but gridded to 512"^. Simulations CI, C2, and 
C3 set Cceii = 1- Simulation C4 uses the 256^ grid with equation 
(HI for Cceii- The additional dumpiness in C2-C4 over CI adds 
structure to the ionization front. Simulation C4 has at least 10 X 
more recombinations than in the other runs. 



in the top left panel (a snapshot from simulation CI, which 
uses a homogeneous density field) has less structure on the 
bubble edges than in the other runs. The picture seen in 
the top left panel is the most similar of all the panels to 
the picture of reionizatio n seen in Monte-Carlo re alizations 
of HII regions usin g the Furlanetto et all (|2004bl l analytic 
model (see figures in lZahn et al. This model spher- 

ically averages the density field around a cell to determine 
its ionization fraction, washing out much of the small-scale 
structure in the density field. The 256"^ and 512"^ runs have 
a similar morphology despite the 512"^ run's higher resolu- 
tion and larger volume- averaged clumping factor. When we 
boost the subgrid clumping factor substantially for the C4 
run, this action does not significantly change the morphol- 
ogy, even though this simulation has a factor of 10 more 
recombinations than the C2 and C3 simulations. 

The A'^x for the Cl-4 runs at Xi^v ~ 0.5 agree to better 
than 20% at scales larger than a Mpc. Simulation CI has the 
least power of all the runs at megaparsec scales because it is 
missing the density-induced structure at the bubble edges. 
In conclusion, the differences in power from clumping in the 
considered models are minor compared to the differences 
that arise owing to the different source prescriptions. 

In the C4 run, the subgrid clumping factor decreases as 
a function of the cell's density (eqn.|4]). We know that over- 
dense regions form substantially more structure, and, there- 
fore, it is possible that the subgrid clumping factor actually 
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increases with density. To test whether such a prescription 
for dumpiness alters the morphology of reionization, we ran 
a small-scale clumping run C5 with Cceii = 4+3 fc, in which 
6c is the baryonic overdensity smoothed at the cell scale. 
This clumping prescription y ields a similar scaling with den- 
sity to the Cceii Pceif that ' Kohler et all (|2005l ) finds in a 
4 Mpc simulation in which the halo particles are also re- 
moved from the density field. This parametrization results 
in a photon to ionized baryon ratio of ~ 2 at the end of 
reionization and Cv ~ 20 throughout reionization. We do 
not plot the results for C5, but we find that the HII regions 
have slightly more structure on the edges in this case than 
in C3 and C4. Overall, the structure of reionization is not 
significantly altered in C5 from the other clumping runs. 

Why does clumping not affect the large-scale morphol- 
ogy of reionization? Qualitatively, large-scale density fluc- 
tuations significantly enhance the mass in sources that are 
present within an overdense region relative to the mean. 
However, the number of absorptions and recombinations per 
unit volume are not enhanced by the same margin. This 
leads to the enhanced abundance of ionizing photons win- 
ning in overdense regions and shaping the morphology of 
reionization. For a more quantitative treatment, one can 
solve for the overdensity that a region must have to be ion- 
ized given some source prescription and parametrization of 
the gas clumping. This overdensity threshold can then be 
used to calculate the bubble size distribution with the ex- 
cursi on set formalism (jPurlanetto et al.ll2004bl : iBond etajj 
Il99ll ) . For reasonable parameterizations of the clumping fac- 
tor, this exercise shows that clumping does not significantl y 
change the bubble morphology for fixed Xi (McQuinn"2006') . 

On smaller scales, density fiuct nations become more im- 
portant in shaping reionization. For a single HII region ion- 
izing a region of 10 Mpc in radius at 2; = 7, the HII region 
is not a perfect sphere, but has fiuct nations in radius with 
AR/R ^ 0.2. These fiuct nations are generated by column 
density fiuctuations between different lines from the source 
to the bubble edges. Lines with lower column densities will 
lead to fingers protruding from the HII regions. Such fea- 
tures are also present when many sources are within a bub- 
ble. 

In addition to imprinting structure on the bubble edges, 
dumpiness has a considerable effect on the part in 10^ fiuc- 
tuations in the neutral fraction within the bubbles. We will 
come back to this in future work. 

In conclusion, quasi-linear density fiuctuations imprint 
substructure on the bubble edges, but do not affect the 
large-scale morphology of the bubbles. Quasi-linear fiuctu- 
ations also increase the number of recombinations and can 
extend reionization. We address the effect of self-shielding, 
non-linear density enhancements in ^6.21 

6.2 Minihalos 

The minimum mass minihalo that retains gas depends on 
the thermal history of the IGM. The Jeans mass at 2; = 10 
for gas that cools adiabatically since thermal decoupling 
from the CMB is 6 x 10^ M© (Barkana & Loeb 2002) 
and the filtering ma ss is approximately ten times larger 
(jGnedin Huil [19981 ). However, reheating by X-rays prior 
to reionization will make the gas warmer than this, eras- 
ing gas density fiuctuations at progressively larger scales. 



iFurlanettol (|2006l l estimates that if POPII stars are respon- 
sible for reionization then the gas temperature is a couple 
hundred degrees Kelvin prior to the time the Universe is 10% 
ionized. This estimate is based on extrapolating local X-ray 
luminosities to high redshifts. A heated, neutral IGM has a 
Jeans mass of Mj =4x10^ M© [T/(200 K) x (l + z)/(10)]^/^ 
An isolated minihalo that holds onto its gas during re- 
heating will subsequently lose its gas via photo-evapo ration 
as ionizing fiux im pinges upon it (Barkana & Lo^ Il999l : 
I Shapiro et al.l 120041' ) . The timescale for photo-evaporation 
tev of a minihalo is roughly the sou nd-crossing time of the 
halo, which for 10^ K gas ionized is ([Shapiro et al. I I2OO4I ) 

This formula works well when the incident fiux is large, but 
under-predicts the evaporation tim e for the ionizing fiuxes 
that are typical during reionization (|lliev et al.ll20059 ). The 
duration of reionization in our simulations is a few hundred 
million years, comparable to the evaporation timescale of 
minihalos (eqn.[5]), suggesting that minihalos will be present 
for all times during reionization. 

Prior to evaporation, a minihalo is optically thick for 
a typical ionizing photon. An incident photon ionizes a hy- 
drogen atom within the minihalo and the photon's energy is 
converted primarily into kinetic energy of the minihalo gas 
rather than into additional IGM ionizing photons. The mean 
free path at z = 6 to intersect a halo of mass (10^, 10^, 10^) 
Mq within a virial radius is (4, 7, 17) Mpc [or at 2; = 12 
is (6, 19, 74) Mpc] if we assume the Press- Schechter mass 
function. 

Several previous calculations have attempted to encap- 
sulate the effect of minihalos via a clumping factor (e.g., 
iHaiman et al.|[2000[ ). We emphasize that this is not an ap- 
propriate way to treat minihalos. Minihalos are self-shielded 
such that the densest inner region s should not contribute to 
the clumping (|lliev et al ] l2005bl ). In addition, in the con- 
text of large-scale simulations, only a small po rtion of pho- 
tons through a cell will intersect a minihalo. ICiardi et al.l 
(2006) was the only previous study to investigate minihalos 
in the co ntext of large - scale radiative transfer simulations. 
However, ICiardi et all (|2006l l set the cell optical depth in 
minihalos to be the average optical depth for all sight-lines 
through the cell. The average optical depth from minihalos 
can be large even though the vast majority of sight-lines will 
not intersect a minihalo. A more appropriate model for the 
minihalos is to treat them as dense absorbers with an ab- 
sorbing cross section amh- We adopt this treatment for the 
minihalos: Only the fraction (Tmh/^ceii of photons in a ray 
that passes through a cell of side-length Lceii are absorbed 
in a minihalo of cross section amh that sits within the cell. 

We add minihalos to our simulation box using the mean 
value method. Method 1 discussed in ^Sl We use the Press- 
Schechter mass function for the minihalos, but using the 
Sheth-Tormen mass function instead would not affect our 
conclusions. The mass function of minihalos is calculated in 
each cell on a 64^ coarse grid, and the mass in mininhalos 
for a coarse cell is divided equally among its fine cells. This 
method is justified because the mean free path for photons 
is always larger than the width of a coarse cell in our models. 

In all of our calculations, we assume that once a re- 
gion is ionized, no new minihalos form owing to "Jeans 
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mass suppression". To incorporate this suppression, we cal- 
culate the opacity of a cell at redshift z that was ionized 
at Zrei from the mass nps{m, 6o,m , Mc, Zrei) rather than 
nps{Tn,So,M, Mc, z). However, we find that our results are 
unchanged if we omit suppression. This is because miniha- 
los are abundant at the redshifts relevant to our study such 
that the number density of minihalos is not rapidly chang- 
ing. For higher redshift reionization scenarios, the degree 
to which minihalos are suppre ssed from formin g in ionized 
regions can play a larger role (ICiardi et al.ll2006l V 

To understand the impact of minihalos, we adopt three 
simplified models for these objects. In our most extreme 
model for minihalos (simulation M3) , we make all minihalos 
with mass greater than 10^ Mq opaque to ionizing photons 
out to a virial radius. The mass cutoff of 10^ Mq IS consis- 
tent with a minimal amount of reheating. Simulation M2 is 
the same as M3, except that the effective cross section (Jmh 
of a minihalo to ionizing photons is not fixed as a function 
of time, but instead the function used for (Jmh is motivated 
by the ev olution of the c r oss se ction in the simulations pre- 
sented in 'Shapir o et al.l (HoOJ) - initially the outer layers 
of the gas in minihalos are quickly expelled leaving a dense 
core, which is evaporated over a time tev The formulas we 
use in M2 for (Jmh and tev are presented in along with 
a discussion of potential drawbacks. Finally, simulation Ml 
has the same sources as the other minihalo runs but does 
not include any minihalos 

Figure [9] plots the ionization history of simulations Ml 
(solid curve), M2 (dotted curve) and M3 (dash- dotted curve). 
All of these simulations use the source luminosity of N(m) — 
9 X 10^^ m/(lO^M0) photons s"^ The absorptions in the 
minihalos extend reionization by less than 100 million years 
in M2 and by more than 250 million years in simulation 
M3. In addition, one in every two ionizing photons in M2 
is destroyed in a minihalo by Xiy — 0.8, and two in every 
three are destroyed in M3. 

Figure [To] shows slices through the Ml, M2, and M3 
simulations (top, middle, and bottom panels, respectively) at 
Xi^v = 0.55 (left panels) and at Xi^v = 0.8 (right panels). 
[Note that, due to a limited number of outputs at which to 
compare, the output for simulation Ml is ^ 7% less ionized 
than the outputs for the other runs.] The total number of 
absorptions inside minihalos increases from simulation Ml 
to M2 to M3. The major effect from minihalo absorptions 
is that the largest bubbles (bubbles larger than the photon 
mean free path) grow more slowly, whereas the growth of 
the smaller bubbles is uninhibited. This effect is particularly 
noticeable in simulation M3, in which the average mean free 
path is 4 Mpc. The mean free path becomes larger than 
this as the smallest halos are evaporated in simulation M2, 
such that the effect of minihalos on the bubble sizes is less 
significant. The smaller bubbles are still larger in M2 than 



8 iBarkana Loebl (l2002l ) finds that minihalos impose a much 
shorter mean free p ath than in our models. The reason for this 
difference is because ISarkana Loe uses a static model 

for the minihalos, which results in each minihalo having a much 
larger cross section. IShapiro et al.l (|2004[ ) finds that the outskirts 
of the minihalo are quickly phot o-evaporated, leaving a smaller 
cross section than in lBarkana Loeb (2002). The parameteriza- 
tions in this section assume the outskirts are quickly evaporated. 
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Figure 9. The volume- averaged ionization fraction for simula- 
tions Ml (solid curve)., M2 (dotted curve), and M3 (dash-dotted 
curve). In M3 the minihalos absorb more photons than in M2, 
and there are no minihalos in Ml. All simulations have the same 
source prescription. The presence of minihalos extends the dura- 
tion of reionization. 



in Ml. (Since Ml is at a ~ 7% smaller x^, if we compared 
at the same x^, this trend would be more noticeable.) 

Figure [TT] shows the bubble PDF for the minihalo runs, 
in which the bubble radius is defined as in 21 We confirm 
that the bubbles are smaller when the minihalos are present, 
particularly once the biggest bubbles become larger than 
the photon mean free path. At xiy — 0.8, the characteristic 
bubble radius is 20 Mpc in Ml (solid curve in Fig. fTT]) . 7 Mpc 
in M2 (dotted curve) , and 4 Mpc in M3 (dot-dashed curve). 
In the minihalo models, the characteristic scale is set roughly 
by the average photon mean free path, which is 4 Mpc in 
simulation Ml. This decrease of the characteristic bubble 
scale fro m the d ense absorbers was first predicted in analytic 
models (jFurlan etto & Oh 2005). However, we do not find the 
sharp cutoff in effective bubble size at the scale of the mean 
free path found in the analytic work of iFurlanetto k, OhI 
(2005). The reasons for this difference are primarily that 
analytic models make the simplifying assumptions that the 
mean free path is spatially uniform and that photons from a 
source cannot travel a distance further than one mean free 
path. 

Figure [12] plots for the Ml (solid curves), M2 

(dotted curves), and M3 (dot-dashed curves) simulations for 
Xi^v ~ 0.55 (top panel) and Xi^v ~ 0.8 (bottom panel). 
The minihalos suppress the large-scale ionization fluctua- 
tions and increase the size of the fluctuations at smaller 
scales. The significance of the effect of minihalo absorptions 
increases with ionization fraction as the bubbles become 
larger. Notice that the total power is contained within the 
box for the models with minihalos in Figure [12] (the power 
peaks at smaller scales than the box scale) - the presence of 
minihalos reduce the size of the box necessary to simulate 
reionization. Note that the differences in A^^, among the 
minihalo models we consider (simulations Ml-3) are not as 
large as the differences in A^^, among the source models for 
Xi = 0.55 (simulations S1-S4, Fig. [5]). However, for larger 
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Figure 10. The effect of minihalos on the ionization maps for 
Xi^Y ~ 0.55 (left panels) and Xi^v ~ 0-8 (right panels). The top 
panels are shces from the Ml simulation in which minihalos do not 
affect the propagation of the ionization fronts (because we have a 
limited number of outputs, the top panels are at approximately a 
7% lower ionization fraction than the others). White regions are 
ionized and black are neutral. The middle panels are from M2, in 
which minihal os are evapora t ed wi th a prescription motivated by 
the results of IShapiro et al.l (|2004l ) and llliev et~aLl (|2005bl ). The 
bottom panels are from simulation M3 in which minihalos are not 
evaporated, and all halos above 10^ Mq absorb ionizing photons 
with impact parameter less than one virial radius (yielding an 
average photon mean free path of 4 Mpc). Minihalos inhibit the 
largest bubbles from growing. 

ionization fractions (see bottom panel) the effect of miniha- 
los on the structure of reionization can be comparable to the 
effect of different source models. 

Dense systems other than minihalos may limit t he mean 
free path of ionizing photons during reionization. iGnedinI 
([2000a) finds that such systems play an important role in 
reionization in radiative-hydrodynamics simulations. The ef- 
fect of these "Lyman-limit" systems should be similar to the 
effect we find for the minihalos. 



7 REDSHIFT DEPENDENCE 

Up to this point, we have only considered reionization 
scenarios in which overlap occurs at 2; ^ 7 and result 



0.5 




1 10 
R (Mpc/h) 



Figure 11. The bubble size distribution for the Ml (solid curve), 
M2 (dotted curve), and M3 (dot-dashed curve) simulations for 
Xi ~ 0.8. For the minihalo models, the bubble size peaks at 
roughly the photon mean free path. 




Figure 12. The ionization fraction power spectrum Aj^, for 
the Ml (solid curves), M2 (dotted curves), and M3 (dot-dashed 
curves) simulations for Xi ~ 0.55 (top panel) and Xi ~ 0.8 (bot- 
tom panel). (Note that the snapshot from Ml is really at a 7% 
lower ionization fraction than the snapshots from M2 and M3 in 
both panels.) At fixed Xi, the minihalos suppress the large-scale 
ionization fluctuations and increase the size of the fluctuations at 
smaller scales. 

in T = 0.06 — 0.08. However, WMAP's measurement of 
r = 0.09 it 0.03 does not rule out overlap at higher red- 
shifts. Further, the popular conclusion that quasar absorp- 
tion spectra require t hat reioniz ation is endin g at z ^ 6.5 
is being hotly debated (IFan et al. 2006; Mesing er Haiman 
'2OO4I: IWvithe Loeb"2 004l : iLidz et al.. .,2006a : IBecker et al. 
2006:l Lidz et a l. 2006b). At higher redshifts, there are fewer 
galaxies above mcooi, enhancing Poisson fluctuations, and 
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Figure 13. The dependence of the morphology on the redshift of 
reionization. Shces through snapshots from the SI and Zl simula- 
tions. Simulation Zl has the same source prescription as SI, but 
with five times the ionizing efficiency. The top panels compare SI 
at z = 8.2 (left) and Zl at 2; = 11.1 (right), both with Xi^v ~ 0-3, 
and the bottom panels compare SI at z = 7.3 (left) and Zl at 
z = 10.1 (right), both with Xi^v ~ 0-6. The ionization fields from 
SI and Zl are strikingly similar. 




Figure 14. Slices through snapshots from the S3 and Z3 simula- 
tions. Simulation Z3 has the same source prescription as S3, but 
with five times the ionizing efficiency. The top panels compare S3 
at 2; = 8.2 (left) and Z3 at 2; = 11.1 (right), both with Xi^y ~ 0.3, 
and the bottom panels compare SI at z = 7.3 (left) and Zl at 
z = 10.1 (right), both with Xi^v ~ 0-6. The ionization field is very 
similar between S3 and Z3. 



the galaxies that do exist are more biased on average. In 
addition, at higher redshifts the Universe is more dense, re- 
sulting in a higher level of recombinations. Finally, at higher 
redshifts the number of galaxies is growing more quickly, 
possibly leading to a shorter duration for the reionization 
epoch. Owing to all these differences, it is interesting to in- 
vestigate how the structure of reionization when comparing 
at fixed Xi changes with redshift. Analytic models predict 
that the bubble size di stribution at fixed Xj is r elatively un- 
changed with redshift (|Furlanetto et al ] l2004bl l 

Figure [T3l compares snapshots from the SI simulation 
and Zl simulation, which has the same sources as SI, but 
where each source has five times the ionizing efficiency. The 
higher efficiency results in reionization occurring earlier by 
a redshift interval of Az ~ 3. The top panels compare SI 
at z = 8.2 (left) with Zl at z = 11.1 (right), both with 
Xi,v ~ 0.3. The bottom panels compare SI at 2; = 7.3 (left) 
with Zl at 2; = 10.1 (right), both with Xi,v ~ 0.6. The 
ionization field is very similar between SI and Zl for fixed 

Xi . 

We also ran simulation Z3, which uses the same source 
prescription as S3 (a oc m?^^), except that the sources in Z3 
are five times as efficient as in S3. More massive sources dom- 
inate the ionizing efficiency in the S3 and Z3 models than 
in SI and Zl. Since the more massive sources are closer to 
the exponential tail of the Press-Schechter mass function, 
the part of the mass function which is rapidly changing, we 
might expect a more significant difference in the ionization 
maps as we change the redshift of overlap than we found in 



the previous case. Figure [TJ] compares the ionization maps 
for the S3 and Z3 simulations (left and right panels, respec- 
tively). The ionization maps are, as with SI and Zl, very 
similar. The differences between the A^^, calculated from SI 
and Zl (or from S3 and Z3) are < 10% at fixed Xi. 

We can understand why the maps look so similar at 
fixed Xi by again comparing the power spectra of the sources 
at these redshifts. The top panel in Figure [2] shows the 
luminosity-weighted power spectrum A^^ of the sources 
used the Sl/Zl simulations (solid curves) and S3/Z3 sim- 
ulations (dotted curves) at z = 6.6 (thick curves) and 11.1 
(thin curves). The differences between A^^ for the SI (or 
S3) sources at z — 6.6 and at z — 11.1 are much smaller 
than the differences between the A^^ for the SI, S3 and S4 
source models. Therefore, we would expect the differences 
between the ionization fields at fixed Xi but separated by 
/\z ~ 3 to be smaller than the differences between the fields 
for the SI, S3, and S4 models, which is what we find. 

Because the ionization maps do not depend strongly on 
the redshift of reionization, we expect that our conclusions 
in previous sections hold for slightly higher redshift reion- 
ization scenarios. The invariance of the ionization fields with 
redshift also implies that the conclusions in this paper are 
not sensitive to the value of as . If reionization occurs at very 
high redshifts, redshifts where the cooling mass sources are 
extremely rare, then the topology of reionization will shift 
from the topology seen in SI to something closer to what 
is seen in S4 - the bubbl es will become larg er and more 
spherical (see discussion in IZahn et al ] (|2006bl )). 
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8 OBSERVATIONAL IMPLICATIONS 

In this section, we briefly discuss the potential of observa- 
tions to distinguish different reionization models. We limit 
the discussion to Lya emitter surveys and 21cm emission. In 
future work, we will discuss the implications for these and 
other observations in more detail. 



8.1 Lya Emitter Surveys 

Narrow band Lya emitter surveys are currently probing red- 
shifts as high as z = 6.5, and projects are in the works 
to search for higher redshifts Lyg emit t ers ( K ashikawa 
I2OO6I : ISantos et al]l2004l : iBarton et aLll2004l : live et al.,.200a ;i . 
If there are pockets of neutral gas at these redshifts, 
the statistics o f these emit t ers can be d ramatically al- 
tered (Madau fc~Reesl I2OOOI: iHaimanl I2OO2I: ISantos 2004 



iFurlanetto et al.ll2006al : iMalhotra Rhoadsll2006h . Sources 
must be in large HII regions for the Lya photons to be able 
to redshift far enough out of the line center to escape ab- 
sorption. Therefore, the structure of the HII regions will 
modulate the observed properties of the emitters. Because 
of this modulation, Lya emitters could be a sensitive probe 
of the HII bubbles during reionization. From the current 
datat on these emitters, there is disagreement as to whether 
there is evidence for reionizati on oX z — 6.5 (jKashikawal 
I2OO6I : iMalhotra Rhoadslliooel ) . 

The calculations in this section are all at z — 6.5, the 
highest redshift at which there are more than a handful of 
confirmed Lya emitters]^ Rather than re-run our simula- 
tions to generate maps with different ionization fractions at 
z = 6.5, we instead use the property that the structure of 
HII regions at fixed xi is relatively independent of the red- 
shift (as demonstrated in 33)- We take the ionization field 
from the simulation for higher z and use this field in combi- 
nation with the z — 6.5 sources. Since the photo-ionization 
state of the gas within an HII bubble is dependent on the 
redshift, we remove the residual neutral fraction within each 
HII region when calculating the optical depth to absorption 
T"Lya • The rcsldual neutral gas primarily affects the blue side 
of the Une, which we assume is fully absorbed. We also ig- 
nore the peculiar velocity field in this analysis. The peculiar 
velocities are typically much smaller than the relative veloc- 
ities due to Hubble expansion between the emitter and its 
HII front, and, therefore, this omission does not affect our 
results. 

Next, we integrate the opacity along a ray perpendic- 
ular to the front of the box from each source to calcu- 
late TLya. Rather than assume an intrinsic Lya line pro- 
file and follow many frequencies, we calculate the optical 
depth TLya for a photon that starts off in the frame of the 
emitter at the line center and set the observed lumi- 
nosity Lobs, Lya = « Lint, Lya exp[-rLya (z^o)] , in which tt is 
a constant of proportionality that encodes the amount of 
absorption at the line center. For reference, an isolated bub- 
ble of 1 proper Mpc that is fully ionized in the interior has 



^ The redshifts that can be probed from the ground are limited by 
sky lines, which contaminate a significant portion of the relevant 
spectrum. At z = 6.5 there is a gap in these lines that allows for 
observations. 



T"Lya(z^o) = 1- Wc also assumc the escape fraction is indepen- 
dent of halo mass such that Lint, Lya = &iVuvE3 In future 
work, we will do a more thorough analysis that includes the 
velocity field, the neutral fraction within the bubbles, as 
well as several frequencies around z/q- We also igno re here 
any stochasticity in the Lya emission from galaxies. ISantosI 
(|2004) discusses the importance of many of the effects that 
are ignored in the calculations in this section. 

Figure [TSlplots the number density of Lya emitters with 
luminosity above aLint,Lya for several volume- averaged ion- 
ization fractions denoted by Xi in the plot. We use the fact 
that there is monotonic relationship between luminosity and 
mass in our models, which allows us to plot mass on the ab- 
scissa. The top panel is from SI in which Lint, Lya oc m and 
the bottom is from S3 in which Lint, Lya cx: m^^^ . Because the 
ionizing sources in S3 are rarer, the bubbles are larger and 
the luminosity function is less suppressed. For both sim- 
ulations, once the Universe is more than half ionized, the 
luminosity function is not significantly suppressed at fixed 
Xi . The normalization of the luminosity function is very sen- 
sitive to ionization fractions Xi < 0.5 in both models. 

The luminosity function for different ionized fractions in 
our calculations is suppressed from the intrinsic luminosity 
function by a factor that is fairly independent of halo mass 
fFig. I15|) . This prediction for the observed luminosity func- 
tion is similar to the analytic predictions of IFurlanetto et aP 
(2006a), which use a similar source prescription to that of 
SI. However, the luminosity function we predict is less sup- 
pressed by a factor of 1.5 — 2. Thi s small difference is partly 
because IFurlanetto et al. 1 (|2006al ) underestimates the free 
path a photon will take inside a bubble. In 'Furlane tto et"al] 
(2006a), for computational convenience the distance for a 
photon to travel within a bubble is defined as the distance 
from the source to the nearest neutral clump rather than 
the distance along the ray to the bubble edge. 

iKashikawal (|200 6) finds significant evolution in the lu- 
minosity function between z = 5.7 and ^ = 6.5 and suggests 
that this might be evidence for reionization. However, the 
z = 6.5 luminosity function differs most with the z = 5.7 at 
the high mass end, as opposed to our prediction of it being 
uniformly suppressed. We suggest that the observed evolu- 
tion is more consistent with cosmological evolution in the 



The precise value of the proportionality constants a and b does 
not matter for the subsequent discussion in this section. The value 
of a and b does matter if we are to compare our results with ob- 
servations. The standard assumption is that a = 0.5 (the blue 
side of the line is absorbed while the re d side is unaf fected), but 
a is probably smaller than this value (ISantosll200l) . In princi- 
ple, we could calculate a from the ionization field in our simu- 
lation, but we leave th is to future work. In the absence of dust, 
b = 0.67 (1 - fesc)huoc (|Osterbrocklll989l ) such that if we assume 
/esc <C 1 then the observed Lya luminosity of these sources is 
La = 3 X 10^^771/(10^° Mq) erg s^^ in simulation SI for a = 0.1. 
The observed emitters have luminosities of 2 x 10^^ — 1 x 10^^ 
erg s"-*^, which correspond to halos with 777 > lO"*^-*^ Mq in SI. 
Presently, surveys cover ~ 10^ Mpc"^ dX z = 6 .5, but probe onl y 
the ~ 100 brightest emitters in that volume (IKashikawal l2006h . 
Assuming for simplicity that all halos host an emitter (which is 
certainly not true in detail), we reproduce th e observed abun- 
dance of Lya emitters n ^ 2x 10"^ Mpc~^ (|Kashikawal l2006h 
if all halos with masses > 3 x lO-*^-*^ Mq host observed emitters 
(assuming as = 0.9). 
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Figure 15. The number density of Lyo; emitters above a certain 
Lya luminosity Lm = aLint,Lya at z = 6.5 and for several dif- 
ferent volume averaged ionization fractions (denoted by Xi in the 
plot). Here we use the fact that there is monotonic relationship 
between luminosity and mass in our source models, which allows 
us to plot halo mass on the abscissa. The top panel is from the 
SI simulation in which I/int,Lya oc m and the bottom is from S3 
in which Lint, Lyo; oc m^/^. Because the sources in S3 are rarer, 
the bubbles are larger and, therefore, the luminosity function is 
less suppressed. Current surveys z = 6.5 probe a volume of 
10^ Mpc^ and have found ~ 100 emitters. 



abundance of massive host halos, rather than reflecting an 
evolving ionized fraction. 

Figure [16] shows maps of the Lya emitters in simulation 
S2 with m > 5 X 10^° M©. This mock survey would subtend 
0.6 degrees on the sky and has a volume of 3 x 10^ Mpc^. 
The left panels are for Xiy — 0.35 and the right are for 
Xi^v = 0.7. The top panels show the average ionization frac- 
tion for a projection of width 31 Mpc, corresponding to a 
narrow band filter with width A A = 100 angstroms. White 
regions are fully ionized and black are fully neutral. The 
middle panels show the intrinsic population of Lya emit- 
ters. There are 1800 of these halos in the survey; the den- 
sity of these halos is an order of magnitude higher than the 
density currently probed by narrow band Lya surveys. The 
bottom panels show the observed emitters [with observed 
luminosity greater than Lint,Lya(m = 5 x 10"^° ^o)], which 
is modulated by the ionization field in the top panel. In 
the left, bottom panel, there are 500 visible emitters and 
in the right, bottom there are 1400. Detecting these large- 
scale variations in the abundance of Lya emitters would be 
a unique signature of patchy reionization. In future work, we 
calculate several clustering statistics from our emitter maps. 

Current surveys at z = 6.5 are dominated by Poisson 
fluctuations and are not yet sensitive to density fluctua- 
tions or bubble-induced fluctuations from reionization. Fig- 
ure [16] illustrates, however, that once surveys resolve enough 
sources then they will be able to detect fluctuations from 
the HII regions (these fluctuations are generally much larger 
than the density- induced fluctuations). At larger scales, the 




Figure 16. Mock z = 6.5 Lyo; emitter surveys at two different 
stages of reionization. The left panels are for Xi^v = 0.35 and the 
right are for Xi^v = 0.7. The panels would subtend 0.6 degrees 
on the sky. All panels use the SI simulation. Top Panels: A map 
of the average ionization fraction from a slice through the box of 
width 31 Mpc. In white regions the projection is fully ionized and 
in black it is fully neutral. Middle Panels: The intrinsic popula- 
tion of Lyo; emitters in the same 31 Mpc slice. Our mock survey 
is sensitive to emitters with halo masses m > 5 x lO-*^*-* Mq. There 
are 1800 such halos in the survey, resulting in a density that is 
an order of magnitude higher than the density probed by current 
Lyo; surveys. Bottom Panels: The observed distribution of Lya 
emitters if the Universe is ionized as given in the corresponding 
top panels. The observed distribution of these emitters is modu- 
lated by the presence of HII regions (see top panels for location 
of HII regions and contrast with intrinsic distribution). 



bubble fluctuations start to dominate over the Poissonian 
fluctuations. The fluctuations generated by the bubbles are 
order unity at the bubble scale, whereas the intrinsic Poisson 
fluctuations are given by A(A;)^ = /c^/(27r^n). If we equate 
this with unity at /c = 0.5 /i Mpc~^ (corresponding to a 
bubble scale of i? = 18 Mpc), we find that surveys must 
be sensitive to source densities of n = 2 x 10~^ Mpc~^ in 
order to overcome Poisson noise and be able to image the 
bubble-induced fluctuations. The requirements for a statisti- 
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cal detection are less stringent. Currently, surveys can probe 
to densities of n ^ 2 x 10~^ Mpc~^. In future work, we pro- 
vide a more quantitative estimate for the number density 
that surveys must probe to detect reionization. 

In future work we will also include the effect of 
minihalos and gas dumpiness on the Lya emitters. 
Minihalos/Lyman- limit systems limit the bubble size and so 
could potentially suppress the observed luminosity function 
more substantially than we find in the SI and S3 simula- 
tions. 



8.2 21cm Emission 

The LOFAR and MWA radio interferometers are being built 
to observe high redshift neutral hydrogen via the 21cm line, 
and the GMRT interferometer can already observe at these 
wavelengths. These telescopes hope to observe an increase in 
brightness temperature over that of the CMB at wavelengths 
A = 21cm (1 -\- z) for z > Zrei with amplitude 
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where Tg is the spin temperature and Sb is the baryonic 
overdensity. Equation (|6]) (as well as our calculations) ne- 
glects redshift-space distortions, which can enhance the sig- 
nal (Barkana & Loeb 2005). However, these distortions offer 
only a small enhancement of the signal on the large scales of 
interest at which ionization fluctuations dominate the signal 
(|McQuinn et al.H 2006). We also assume Ts ^ Tcmb in this 
section, likely a good approximation during the bulk of the 
reionization epoch (Ciardi & Madau 2003; Furlanetto 2006*). 

Figure \T7\ plots the 21cm power spectrum for the SI 
(solid lines), S2 (dashed lines), S3 (dot-dashed lines), and S4 
(dotted lines) simulations for Xi^v = 0.2 (top panel), xiy — 
0.5 (middle panel), and Xi^v — 0.7 (bottom panel). The S3 
and S4 simulations have much more power at large scales 
than the other runs, particularly at early times owing to the 
larger bubbles in these runs (Fig. 2]). The signal is very flat 
on the scales probed by the box for most xi. If we had a 
larger box, a sharp decline in power would be observed at 
larger scales than the bubbles. 

The power spectra in Figure [TTl do not include absorp- 
tions from minihalos. If minihalos are as abundant in reality 
as they are in models M2 and M3, this will suppress sig- 
nificantly the large-scale power in the 21cm signal (see the 
in Fig. I12p . The effect of the minihalos on the 21cm 
power spectrum is qualitatively different from the effect of 
changing the sources and should also be observable. 

The projected 1-cr errors for MWA for a 1000 hour ob- 
servation in a 6 MHz band in bins of width A/c = 0.5 /c 
are shown in the middle panel in Figure 1171 The sensitivity 
of LOFAR is comparable to that of MWA. The details of 
this se nsitivity calculation are discussed in iMcQuinn et al.l 
(|2006l V Because of foregrounds, experiments will encounter 
difficultly detecting s maller /c-modes than are plotted here 
(|McQ uinn et al ][20Q6). The first generation of interferome- 
ters are most sensitive to k greater than 0.1 /i Mpc~^ and 
less than 1 h Mpc~^ MWA and LOFAR should be able to 
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Figure 17. The 21cm power spectrum for the SI (solid curves), 
S2 (dashed curves), S3 (dot- dashed curves), and S4 (dotted 
curves) simulations [A2i(k) = A:^(T2i (/c)^}/27r^]. For the top pan- 
els, Xi^v ~ 0-2 (xi^M ~ 0-3), for the middle panels, Xiy ~ 0.5 
(xi,M ~ 0-6), and for the bottom panels, Xi^y ~ 0.7 (xi^M ~ 0-8)- 
In S2 the lowest mass sources dominate (with masses m ~ ''TZcool)? 
and in S4 the highest mass sources dominate (m ~ 5 x IQ-'^'-' Mq). 
At scales k <lh Mpc"-*^, A'^-^ scales approximately as Aj^, such 
that the 21cm signal is a sensitive probe of the bubble structure. 
The error bars are the detector noise plus cosmic variance errors 
on the power spectrum for MWA, assuming 1000 hours of inte- 
gration and a bandwidth of 6 MHz. Foregrounds will eliminate 
the sensitivity to the signal for /c < 0.1 Mpc"-*^. 



distinguish between the SI, S3 and S4 reionization scenar- 
ios at a fixed ionized fraction. If we marginalize over the 
ionized fraction, it is unclear whether MWA can still distin- 
guish between these models. The second generation 21cm 
experiments MWA5000 and SKA will be at least 10 x m ore 
sensitive than MWA and LOFAR (|McOuinn et al.ll2006h . 

In addition, it might be possible to use the evolution 
of the 21cm signal to separate models. For the models con- 
sidered in this paper, the duration of reionization is fairly 
short, spanning an interval of = 2 — 4. It is quite possible 
that upcoming 21cm experiments will be able to observe the 
entire breadth of this epoch. The duration of reionization is 
shortest if the largest mass sources dominate the ionizing 
budget. Also, minihalos tend to cause a delay at the end of 
reionization ( Fig.[9|). Perhaps combining information on the 
duration of reionization with the power spectrum at different 
times can help us understand the source properties as well 
as the role of the minihalos. Higher order terms in the 21cm 
power spectrum may aid in distin g uishing different re ion- 
ization scenarios (|Lidz et al ] |2006d ). IZahn et al.1 (|2006al ) (in 
preparation) investigates how well upcoming 21cm experi- 
ments can constrain certain reionization models. 



20 M. McQumn et al 



9 CONCLUSIONS 

We have run a suite of 94 Mpc radiative transfer simula- 
tions to understand the size distribution and morphology 
of HII regions for 0.1 < Xi < 0.8. These simulations are 
the first that include sources down to the cooling mass and 
that are large enough to contain many HII regions. We have 
incorporated structures that all large-scale simulations of 
reionization do not resolve with analytic prescriptions. 

We find that the morphology of HII regions is most 
sensitive to the parameter Xi. If we compare different reion- 
ization scenarios at the same x^, they tend to look similar. 
This is not to say other factors besides Xi do not change the 
morphology. The sources responsible for reionization are the 
second most important factor. If we compare at fixed Xi, we 
find that the HII regions become larger (by as much as a 
factor of 4) and more spherical as the sources become rarer. 
The bubbles are larger for the rarest sources because these 
sources are the most biased. 

The next most important factor for shaping the mor- 
phology is the presence of minihalos. Once the mean free 
path for a photon to intersect a minihalo becomes smaller 
than the bubble size, the effect of minihalo absorptions be- 
comes important. As a result, minihalos inhi bit the largest 
bubbl es fromgrowin^ I f we us e the results of IShapiro et al.l 
(|2004l l and llliev et aL 1 (|2005b) to characterize the miniha- 
los, we find that these objects have a modest effect on the 
overall properties of the HII regions at fixed x^, decreasing 
A XX by as much as 50% for the largest modes in our box. 
In a more extreme case we considered, in which the average 
mean free path is 4 Mpc during reionization, the impact of 
minihalos is even more substantial. Minihalos do not have 
the same effect as changing the source efficiency. 

We find that thermal feedback and quasi-linear den- 
sity inhomogeneities have more minor consequences for the 
topology of the bubbles at fixed Xi . This is fortunate because 
these quantities are poorly constrained. Feedback does not 
substantially change the morphology of reionization at fixed 
Xi because the bias difference between the m > 10^ Mq 
halos and the m > 10^ Mq halos is relatively small. (The 
typical halo that is suppressed by feedback is located in a 
similar region as the typical halo which is not.) Megaparsec- 
scale, quasi-linear density fluctuations add structure to the 
boundaries of the HII regions. This additional structure is 
ignored in analytic models. However, as we increase the level 
of small-scale gas clumping, either by increasing the resolu- 
tion or by increasing the subgrid clumping factor, the large- 
scale structure of the HII regions is largely unaffected at 
fixed Xi . This is true even if small-scale gas clumping results 
in a substantial number of recombinations. We find that the 

at fixed Xi differ by no more than 20% as we vary the 
volume averaged clumping factor from to 30. The qualita- 
tive reason why clumping does not affect the morphology of 
reionization at fixed Xi is because the enhanced photon pro- 
duction in a large-scale overdense region (that is a bubble) 
is always able to overcome the enhanced number of recom- 
binations, even in extreme clumping models. 

The conclusions in this paper hold if overlap occurs at 
slightly higher redshifts then in our typical simulation in 
which 2;overiap ~ 7. lu fact, we found that if we boosted 
the source efficiencies such that at ^overlap ~ 10, the ion- 
ization map is essentially unaffected. We showed that this 



can be explained by the relatively small differences in the 
luminosity-weighted source power spectrum at z = 7 com- 
pared to that at 2; = 10 in the models we considered. The 
conclusion that the structure of reionization does not de- 
pend on the redshift is no longer true if we compare with 
simulations that reionize at much higher redshifts, redshifts 
at which the sources become extremely rare. In this case, 
reionization may have a similar morphology to simulation 
S4, in which the rarest sources dominate. 

In this paper, we did not concentrate on predicting the 
duration of reionization. However, many of the effects we 
consider impact the duration of reionization, even if they do 
not impact the morphology. We find that our most extreme 
minihalo model extends the duration of reionization by 250 
million years {Az ~ 1.5). In addition, feedback on POPII- 
like ionizing sources from photo-heating can in extreme cases 
extend reionization by 200 million years. 

Analytic models provide a convenient and intuitive 
framework to understa n d th e structure of reio nization 
(iFurlanetto et al.1 l2004bl . l2005l : IZahn et all l2006bh . These 
models do not suffer from the same scale limitations as sim- 
ulations, and they supply a quick method to explore the 
large parameter space relevant to reionization. In addition, 
these models enhance our physical intuition regarding the 
processes that shape this epoch. We have confirmed the an- 
alytic predictions that the bubble size distribution is approx- 
imately log-normal and that the sizes of the bubbles increase 
as the sources become more biased. Further, we confirm the 
prediction of analytic models that bubble sizes are largely 
unchanged if we compare the same model at different red- 
shifts, yet fixed ionized fraction. We also showed, however, 
that current analytic models encounter some difficulties in 
describing the effect of minihalos and of Poisson fluctuations 
in the source abundance on the structure of reionization. An- 
alytic models cannot incorporate the sophisticated models 
for thermal feedback, gas clumping, and minihalo evapora- 
tion that it is possible to include in radiative transfer simu- 
lations. 

Upcoming observations have potential to distinguish 
the source models we considered. We make predictions for 
the luminosity function of Lya emitters as a function of 
Xi. We construct maps of Lya emitters from a mock sur- 
vey that show large-scale fluctuations in the distribution 
of emitters due to the HII regions, suggesting that future 
measurements of the clustering of emitters may reveal the 
signature of patchy reionization. Future 21cm arrays hold 
much promise for probing reionization; measurements of the 
power spectrum with the MWA and the LOFAR arrays can 
distinguish the SI, S3 and S4 source models. 

Upcoming observations can reduce the parameter space 
that reionization simulations need to explore. If we can mea- 
sure the number of ionizing photons produced by high mass 
galaxies and bright quasars at high redshifts, this will reduce 
the almost total freedom we currently have in the ionizing 
luminosity. Observations of the mean free path of ionizing 
photons at high redshifts may reveal whether the Lyman- 
limit systems are the minihalos as well as how fast these 
systems are being evaporated. A precise measurement of 
the Thomson scattering optical depth from the CMB will 
constrain the average redshift of reionization. 

It is important to continue to improve large-scale sim- 
ulations of reionization to understand the reionization pro- 
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cess in more detail. Future simulations need to investigate 
the effect of more realistic star formation models, metal pol- 
lution, and alternative sources of ionizing photons. In addi- 
tion, larger simulations than are presented here are neces- 
sary to statistically describe this epoch for x^ > 0.7. It is 
also useful to run small-scale simulations to develop more 
realistic subgrid parameterizations for the minihalos and for 
the clumping factor. These parameterizations will be essen- 
tial for modeling the end of reionization, a time when the 
rate of evaporation of the Lyman-limit systems plays a key 
role in determining the structure of reionization. Also, such 
parameterizations are necessary to extend our calculations 
to model accurately the part in lO'^ neutral fraction fluctu- 
ations that characterize the high redshift Lya forest. 

An accurate interpretation of future observations of 
reionization, while certainly challenging, does not appear 
impossible. The characteristics of HII regions during reion- 
ization might have depended on a huge number of poorly 
constrained parameters, making it impossible to interpret 
observations of this epoch. We find that this is not the 
case. The morphology of the HII regions at fixed Xi boils 
down primarily to the properties of the sources and of the 
minihalos/Lyman- limit systems. 
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APPENDIX A: RADIATIVE TRANSFER 
ALGORITHM 

For the simu l ations in t his p aper, we employ the 
'Sokasian et ajj (|200lL l2QQ3l . [2004) cosmological radiative 
transfer code, but with several significant changes that are 
discussed below. This algorithm inputs grids of the den- 
sity field as well as a list of sources and then casts rays 
from every source, randomizing the order of the sources 
within this loop. Radiative feedback on the density field is 
ignore d. Rays are split adapt ively using the HealPIX algo- 
rithm (jAbel Wandeltll2002l ) such that, at a minimum, N 
rays from a source intersect every cell face (for this paper, we 
set N — 2.1). This algorithm does not iterate the ray casting 
within each time step to converge to the correct ionized frac- 
tion in each cell. Instead, once a cell has been ionized by a 
source within a time step, rays from other sources will pass 
through it. This simplification allows for the algorithm to 
process more sources and larger volumes than other codes. 
In the limit of few sources and few time steps, this simpli- 
fication can lead to artificial structure in the HII regions. 
However, with the vast number of sources in the simula- 
tions in this paper, even with relatively coarse time steps 



we choose, this artificial structure is minimized (as we will 
demonstrate later in this section). 

The temperature history of the gas is not tracked by 
this code. Instead, the code assumes that ionized regions 
are at T = 10^ Kelvin. The temperature affects the number 
of recombinations in the simulation because as oc T~°"^, as 
well as the detailed photo-ionization state of the gas within 
the HII regions. The analysis we have done in this paper 
does not depend on the photo-ionization state of the gas. In 
addition, the value for the subgrid clumping factor, which 
determines the number of recombinations with in a cell, is 
highly uncertain, such that we would not gain any precision 
from including a full calculation of the IGM temperature. 

What follows is a list of the important modificat ions 
that we have made to the original ISokasian et al.1 (|200lh al- 
gorithm: 

• Previously, cells were either ionized or neutral. Cells can 
now be fractionally ionized. We assume that the ionizing 
front is paper thin such that each cell can be broken up into 
a neutral part and an ionized part. 

• Each ray holds a number of photons. In the original 
I Sokasian et al.l (|200lh algorithm, the first ray that hit a cell 
from a particular source carried all the information that the 
cell needed about the source. Subsequent rays from the same 
source did not affect the cell. The advantage to having each 
ray contain a specific number of photons is that it is trivial 
to conserve photons as well as to include photon sinks. The 
disadvantage is that the ionization front has a numerical 
width that is wider than in the previous algorithm. We find 
that the width of the front in the new algorithm is approxi- 
mately 2 cells for a single source. The thickness of the front 
is smaller than 2 cells in the limit relevant to this paper of 
many faint sources. 

• The orientation of the HealPIX ray casting scheme is 
randomly rotated between each time step, and the order 
with which the rays are initially cast is also randomized. 
When rays split adaptively, the order is again randomized 
over the daughter rays. All of this randomization is done to 
minimize artifacts owing to the order in which operations 
are performed. 

• Once a ray has traveled a distance equal to 
-Rbox, proper [3iVsoiirce/(47riVbox)]"^^^, it cau no longer split 

into daughter rays, where iVsource is the ionizing luminos- 
ity from the source and A^box is the total luminosity of all 
the sources in the box. We set 77 = 5 for this paper. Until 
this distance, rays split adaptively such that a set number of 
rays intersect every cell face. This simplification is justified 
by the large numbers of sources in an HII region (typically 
more than 1000 sources), making it unnecessary to have rays 
from one side of an HII region cover the entire front on the 
other side. Our approximation results in the correct fluxes in 
the cells in the limit of many sources. We have investigated 
quantitatively whether this simplification makes a difference 
in the ionization maps. The middle panel in Figure [A2l plots 
the cross correlation coefficient at two times between a sim- 
ulation with no ray termination and a simulation with the 
prescription for ray termination used in this paper (see the 
caption in Fig. IA2l for the definition of the cross correlation 
coefficient). There is essentially no difference between the 
maps. This simplification results in the algorithm running 
over a factor of 5 faster at high ionized fractions. 
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• The previous algorithm reset the density in each ceh 
after a time step to the density field in the next snapshot, 
but did not change the ionized fraction in the cell to ac- 
count for the dynamics of the ionization field. For example, 
a cell that becomes fully ionized would remain fully ionized 
in subsequent time steps (neglecting recombinations), even 
if it gained neutral material from a neighboring cell during 
these time steps. This resulted in the total number of ion- 
ized atoms not being conserved by the previous algorithm 
between time steps. To remedy this issue in the current al- 
gorithm, we account for a dynamic density field by assigning 
some ionization fraction to each particle in the N-body simu- 
lation and then re-gridding the ionization map between time 
steps to account for the particle dynamics. We suspect that 
other cosmological radiative transfer algorithms performed 
on top of a static density grid ignore the dynamics of the 
ionization field in their computations 

• N-body particles that are associated with halos are not 
included in the density field used by the radiative transfer 
algorithm. Otherwise, cells with sources would have substan- 
tial overdensities, and ionizing photons from within the cell 
would have to ionize these cells prior to escaping into the 
IGM. These absorptions are already counted in the escape 
fraction. Removing the halo particles from the gridded den- 
sity field is also appropriate for rays incident on this cell. 
The gas in the massive source halos has cooled to form a 
small disk that is much smaller than the cross section of the 
cell. Therefore, the vast majority of photons coming from 
exterior to the cell do not intersect this disk. The gas within 
galaxies during reionization absorbs a negligible amount of 
external photons because the mean free path of these pho- 
tons to intersect a galaxy is large (larger than the 94 Mpc 
box size employed in this paper) . 

We subjected the radiative transfer algorithm to sev- 
eral tests. As a simple test, we put one source with N — 

photons s"^ in a 65.6 Mpc/h box with 256^ cells, with 
each cell at the mean density of the z = 6 Universe, and set 
the clumping factor C in each cell to C = 1 or C = 30. In 
Figure IA1[ we compare the fraction of the box that is ion- 
ized in this test to the fraction that is predicted by theory 
(using coarse time steps of 5 x 10^ years). Even with such 
coarse time steps, this algorithm matches the theory curves 
well. 

Because the algorithm does not iterate to find the ion- 
ized fraction, this might lead to artificial structure if the 
time step is too coarse. In the limit of an infinitely small 
time step, this algorithm gives us the exact solution. In this 
paper we use a time step of 50 million years. To test con- 
vergence, we run two cosmological simulations on the 256^ 
grid (which we label simulation 1 and 2), using the halos 
with m > 2 X 10^ M© as our sources. Each simulation uses 
a different set of random numbers to establish the order 
of the sources for ray casting. If the 50 million year time 
step is too coarse, the ionization maps from these two sim- 
ulations would differ substantially, whereas the time step 
is sufficiently small if the ionization maps differ insignifi- 
cantly. The bottom panel of Figure IA2I plots the cross cor- 
relation coefficient r = Px-^x^/ \/ Px-^x-i Px2X2 between these 



Note that our code still ignores thermal feedback and therefore 
does not capture the full dynamics of the gas. 
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Figure Al. The fraction of the simulation box that is ionized by 



a single source with = 10' 
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This calculation is performed 



at 2; = 6 with all cells at the mean density of the Universe and 
with the subgrid clumping factor C equal to 1 and 30. The solid, 
dashed and dotted curves are from a theoretical calculation for no 
recombinations, recombinations with C = 1 and recombinations 
with C = 30, respectively. The pluses and crosses are the radiative 
transfer algorithm with C = 1 and C = 30 with coarse time 
steps of At = 5 X 10^ years. The radiative transfer agrees well 
with the theoretical result, slightly under-predicting the number 
of recombinations. 



two runs for ionization fractions of Xi = 0.2 {solid curves) 
and Xi = 0.7 (dashed curves). The cross correlation coef- 
ficient is close to unity at most scales, dropping to 0.8 at 
the cell scale. Note that the cross correlation coefficient is a 
stringent test. We have also looked at the power spectrum 
of these runs. The power spectrum of the ionized fraction 
differs negligibly between these two runs, differing by about 
0.3% at /c = 10 /i Mpc-^ and 1.5% at /c = 20 /i Mpc-^ 

We have also investigated whether the ionization struc- 
ture seen in our fiducial runs has converged to the true struc- 
ture by either increasing the number of rays that are cast 
or by increasing the mesh size to 512^. First, we ran with 
a simulation that casts a much larger number of rays per 
source than the fiducial number of rays (64 x more initially 
and with the criteria that a minimum of 4.1 rather than 2.1 
rays intersect every cell) . We find that the fiducial number of 
rays is enough to capture the structure of the HII region (top 
panel, Fig. IA2|) . Second, we have run a resolution test, com- 
paring a higher resolution 512^ radiative transfer simulation 
without recombinations to a 256"^ run without recombina- 
tions. We re- grid the 512^ snapshots to 256^ resolution for 
comparison. We find that for all Xi the power spectra have 
converged to within 10% at scales with k < 10 h Mpc~^. 
The agreement is even better then this for Xi > 0.5. 

The computation speed of this code scales in a time step 
as Nr Ns where Nr is the number of rays through each 
cell, R is the characteristic bubble size, Ns is the number of 
sources. When the bubbles are large, the code slows down 
considerably. We are working on ways to ameliorate this 
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Figure A2. The cross correlation coefficient r = 
i^ccicc2(^)/\/^i^7(^)~^2^2W between the ionization 
maps from different simulations. All panels are run with 
N{m) = 3 X 10^9 m/(10^ Mq) and using only halos with 
m > 2 X 10^ Mq. The cell scale is at 20 /i Mpc~^ Top 

Panels: The cross correlation coefficient for a simulation that 
casts 64 times the number of initial rays and ensures that 2 
times more rays intersect every cell as in the fiducial runs. 
The solid curve is for Xi = 0.2 and the dashed is for Xi = 0.5. 
Middle Panels: The cross correlation coefficient of a simulation 
that terminates rays with the condition given in this section 
versus a simulation with no ray termination. The solid curve 
is for Xi = 0.3 and the dashed is for Xi = 0.5. The maps are 
extremely similar even though the simulation without ray casting 
took five times longer to reach Xi = 0.5. Bottom Panels: Cross 
correlation coefficient for two runs with two different sets of 
random numbers. The random numbers set the order of the 
sources for ray casting. The solid curve is for Xi = 0.2 and the 
dashed is for Xi = 0.7. All panels show that the ionization field 
has converged sufficiently in our simulations. 



issue with the code. Most simulations in this paper took less 
than two days to reach = 0.8 on a single 2.2 GHz AMD 
Opteron processor. 



APPENDIX B: FITTING FORMULA FOR 
MINIHALO EVAPORATION 



llliev et aT (2005b) provide fitting formula for the evapora- 
tion of minihalos by POPII stars. These simulations do full 
radiative hydrodynamics on minihalos, which are modeled 
prior to front crossing as truncated isothermal spheres (TIS) 
with self similar infall. They provide the formula for the 
evaporation timescale 



150 



M 



0.1 + 0.9 



10 



^-0.35+0.05 logio(i^) 



Myr, 



(Bl) 



where F is the flux (which is time-independent in their sim- 
ulations). To apply these formula to simulation M2, we use 
for F the time averaged flux incident on a cell, with averag- 
ing starting after t he cell becomes ionize d. 

In Figure 29 in lShapiro et al.l (|2004h , the effective cross 
section of a halo for absorbing an ionizin g photon as afunc- 
tion of time is plotted for a lO'^ M© halo. llliev eraD (|2005bl l 
does not provide parameterized fits to the effective cross sec- 
tion, which we need in our calculations. To pro ceed, we fit by 
eye th e curve for the effective cross section in lShapiro et alj 
(|2004l ). We find the function 



C"mh 



X 10 



' ( tev ) 



(B2) 



where n = 0.754 [M/(10^ Mo)]^/^ 10/(1 + z) is the scale 
radius for the TIS profile. By construction in the simula- 
tions of IShapiro et al. = 7rr? at t = 0. However, 
on a timescale of order a million years the outskirts of the 
halo are evaporated, consuming a meager amount of pho- 
tons. The denser inner regions take a significantly longer 
time to evaporate. We set (Jmh = Trr? for the first 5 mil- 
lion years, and subsequently use equation ()B2|) in run M2. 
Of course, the function we use for amh likely does not scale 
correctly with redshift or halo mass. We anticipate that it 
over-predicts the cross section for halos with m < 10^ M©, 
since the gas in the outskirts of these smaller halos will be 
easier to evaporate. 



APPENDIX C: FILTERING MASS 

If all of the gas in the IGM is cold before reionization and 
subsequently jumps to 10^ K at arei, the filtering mass is 
(jGnedin Huilll998l l 



Mf = Mj[- 



l + 4(^)'''-5(^) 



3/2 



(Ci) 



This mass scale can be much smaller than the Jeans mass 
for 10^ K gas Mj and is typically time dependent. Since 
Mf typically corresponds to non-linear scales where linear 
theory is a poor approximation, it is unclear how well equa- 
tion ()Cip represents the smallest mass at which the baryons 
clump. However, Gnedin & Hui (1998) smooth N-body sim- 
ulations by including a pressure force that becomes impor- 
tant at the filtering scale. They conclude that this procedure 
reproduces well the small-scale gas power spectrum seen in 
hydrodynamics simulations. Furthermore, iGnedinI (|2000bh 
finds that the filtering mass provides a good fit to the min- 
imum formation mass for a gas-rich halo. 
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